Cosmin Burtea
email: cburtea@imj-prg.fr

Sergey Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

Charlotte Perrin
email: charlotte.perrin@univ-amu.fr

Hamilton's principle of stationary action in multiphase flow modeling

These lecture notes are concerned with the derivation of the fluid mechanics equations via Hamilton's principle of stationary action. We recall the main conceptual tools of this variational principle which originally applies to classical finite-degrees-of-freedom mechanics and we explain how these tools can be adapted in a continuous framework, in particular for the derivation of the well-known Euler equations describing the motion of inviscid fluids. The core of these notes is the application of Hamilton's principle to multiphase flows. We present a new Lagrangian point of view for the derivation of two-phase flow equations.

Introduction

Hamilton's principle of stationary action is a variational principle that allows one to obtain the equations of motion (Euler-Lagrange equations) for a given mechanical system. The action is the integral over a finite time interval of the corresponding Lagrangian which is the difference between the kinetic and potential energy of the system. The corresponding Euler-Lagrange equations can be seen as 'Newton's laws' governing the mechanical system. In a nutshell, this principle states that the real trajectory of a system between an initial and final configuration in a specified time is found by selecting from all possible trajectories the one for which the first variation of the action vanishes. In some cases, additional constraints (geometrical and physical) may also be imposed for the class of variations admissible variations. The main advantages of the variational point of view are :

• the whole physics is contained in the definition of a scalar function -the Lagrangian of the system.

• due to the Noether theorem, the fulfillment of the basic physical conservation laws is guaranteed by this approach.

• numerical methods and homogenization techniques can be effectively developed for the corresponding Euler-Lagrange equations which is useful for engineering problems.

The main purpose of these lecture notes is to present some recent developments of such a method used in order to derive equations governing multiphase flows.

Organization The present notes are divided into four parts: in the first part dedicated to standard discrete mechanical systems, we present the formalism of Hamilton's principle yielding the well-known Euler-Lagrange equations; the second part of the notes is devoted to the extension of the previous concepts to the continuous (infinite dimensional) systems. This allows us to recover in the third part the classical (compressible and incompressible) Euler equations for a single fluid phase. Finally, the last section presents the application of Hamilton's principle to the derivation of fluid equations describing the dynamics of multiphase flows.

These notes are intended to be an introduction to multiphase flow modelling via variational principles and are addressed to wide audience: graduate/post graduate students, applied mathematicians, physicists etc.

Before introducing Hamilton's principle, let us recall some notations and formulas that we shall use in the rest of the notes.

Notations and useful formulas

In the computations that follow we will use the Einstein convention of summation over repeated indices.

Matrices For any k, n ∈ N * , we denote by M k×n (R) the space of matrices (or second order tensors) with k lines and n columns considered in an appropriate basis. For the discussion below, let us fix a matrix A = (a ij ) i,j∈ 1,n ∈ M n×n (R). The trace of A is the sum of the diagonal elements of A and we denote it by trace A = a ii .

For any i, j ∈ 1, n , we denote by M ij (A) the determinant of the (n -1) × (n -1)-type matrix obtained by removing the i-th line and the j-th column of A. This quantity is referred to as the (i, j)-minor of A. We recall Laplace's formulae

a i0j (-1) i0+j M i0j = a ij0 (-1) i+j0 M ij0 = det A, a i0j (-1) i1+j M i1j = a ij0 (-1) i+j1 M ij1 = 0. (1.1)
which hold for all i 0 , i 1 , j 0 , j 1 ∈ 1, n with i 0 = i 1 and j 0 = j 1 . The adjugate matrix of A, denoted adj (A) ∈ M n×n (R) is given by (adj (A)) ij = (-1) i+j M ji .

Recall that any A ∈ M n×n (R) with det A = 0 is invertible and Laplace's formulas give us

A -1 = 1 det A adj (A) .
Differential calculus Let n, k ∈ N * and Ω be a smooth (say C 1 ) bounded domain of R n . For any φ : Ω ⊂ R n → R k the differential at x ∈ Ω is the matrix Dφ |x ∈ M k×n (R) given by

Dφ |x ij = ∂ j φ i (x)
for all i ∈ 1, k and j ∈ 1, n . Note that, whenever it is not ambiguous, we shall remove the x from the notation of the differential. We denote the gradient by ∇φ = (Dφ) T .

Let Ω 0 be a smooth bounded domain of R m , m ≥ 1, and

f : Ω 0 → Ω. Then D (φ • f ) ∈ M k×m (R)
and for all x ∈ Ω 0 we have the chain rule

D (φ • f ) (x) = Dφ |f (x) Df |x . (1.2)
From now on, m = n and we assume that there exists

X ∈ C 1 (Ω 0 , R n ) ∩ C 0 (Ω 0 ; R n ), which is a bijection between Ω 0 and Ω such that X -1 ∈ C 1 (Ω; R n ) ∩ C 0 (Ω; R n ).
Of course such an X is also a homeomorphism between the boundaries of the two domains and obviously ∂Ω = X (∂Ω 0 ) .

For any function

φ : Ω → R k we denote by φ : Ω 0 → R k the function φ = φ • X.
(1.3)

Using the chain rule (1.2), we have that

DX -1 |X(x) DX |x = I n ∀ x ∈ Ω 0 ,
which rewrites using the convention (1.3)

(DX) -1 = DX -1 . (1.4) Again, owing to (1.2) we have D φ = Dφ DX,
and consequently Dφ = D φ (DX) -1 .

(1.5)

Remark 1.1 Anticipating Section 3, the diffeomorphism X will represent a fluid particle trajectory and the notation • will correspond to a Lagrangian description of the dynamics, that is a description which is attached to the particle trajectories.

Useful formula In particular, formula (1.5) applied to a vector field

B : Ω → R n yields div B = ∂ i B i = trace DB = trace D B (DX) -1 = ∂ k B i (DX) -1 ki . (1.6)
For a real-valued function b : Ω → R, we are interested in obtaining a formula for ∇b. On the one hand, Equation (1.5) holds. On the other hand, fixing det DX > 0, we have for any test function

ψ ∈ C 1 c (Ω; R n ) Ω ∂ i b (y) ψ i (y) dy = - Ω b (y) ∂ i ψ i (y) dy = - Ω0 b (x) div ψ (x) det DX (x) dx = - Ω0 b (x) ∂ k ψ i (x) (DX) -1 ki det DX (x) dx = Ω0 ∂ k (DX) -1 ki b det DX (x) ψ i (x) dx.
Since, by a simple change of variable,

Ω ∂ i b (y) ψ i (y) dy = Ω0 ∂ i b (x) ψ i (x) det DX (x) dx
we deduce the following "conservative formula"

∂ i b det DX = ∂ k (DX) -1 ki b det DX = div (DX) -1 e i b det DX ∀ i ∈ 1, n , (1.7)
where e i is the i-th vector of the canonical base and where we have defined the divergence of the second order tensor by taking the divergence of each column (sometimes in the literature one uses an opposite definition by taking the divergence of each line).

Identifying M n×n (R) with R n 2 one can see the determinant as a function from R n 2 to R. Laplace's formula allows us to compute its differential. More precisely, using (1.1) we see that for all i 0 , j

0 ∈ 1, n ∂ det ∂a i0j0 (A) = (-1) i0+j0 M i0j0 . Moreover, if A = A(s) ∈ C 1 (R; M n×n (R)
), then we have that

d det A ds (s) = ∂ det ∂a ij (A (s)) ∂a ij ∂s (s) = (-1) i+j M ij (A (s)) ∂a ij ∂s (s) = det A (s) A -1 ji ∂a ij ∂s (s) = (-1) i+j M ji (A (s)) ∂a ji ∂s (s) = adj A (s) ∂A ∂s (s) ii = trace adj A (s) ∂A ∂s (s) .
In particular, if det A (s) = 0 then we may write that

d det A ds (s) = det A (s) trace A -1 (s) ∂A ∂s (s) . (1.8)
2 The classical formulation of the principle of stationary action

Consider a system of N interacting particles which move in a three-dimensional space. We assume that the evolution of the system can be completely characterized by the knowledge of:

• generalized coordinates X (t) = (X 1 (t) , X 2 (t) , • • • X N (t)) ∈ R 3N where, for i ∈ 1, n , X i (t) = X 1 i (t) , X 2 i (t) , X 3 i (t) ∈ R 3
represents the position of the i-th particle at time t.

• generalized velocities Ẋ(t) = Ẋ1 (t), Ẋ2 (t), • • • ẊN (t) ∈ R 3N where, for each i ∈ 1, n , Ẋi (t) = Ẋ1 i (t) , Ẋ2 i (t) , Ẋ3 i (t) ∈ R 3 is the time derivative of X i (t)
and it stands for the velocity of the i-th particle at time t.

The evolution equations of the system will be encoded in a given function L = L(t, X(t), Ẋ(t)) called Lagrangian. We fix two times t 0 < t 1 ∈ R + and define the associated Hamiltonian action between t 0 and t 1 by

A t0,t1 [X] := t1 t0 L t, X (t) , Ẋ (t) dt.
(2.1)

We will denote for each i ∈ 1,

N ∂L ∂X i = ∂L ∂X 1 i , ∂L ∂X 2 i , ∂L ∂X 3 i , ∂L ∂ Ẋi = ∂L ∂ Ẋ1 i , ∂L ∂ Ẋ2 i , ∂L ∂ Ẋ3 i .
Suppose that we know the initial and final configurations of the system, i.e. X (t 0 ) = X 0 and X (t 1 ) = X 1 for some given X 0 , X 1 ∈ R 3N . We say that {X(•; µ)} µ∈(-µ0,µ0) is a family of virtual motions (or virtual displacements) if

• X (t; 0) = X (t) for all t;

• X (t 0 ; µ) = X 0 and X (t 1 ; µ) = X 1 .

Therefore, X(•, µ) can be seen as a perturbation of X that has the same initial and final configurations as X (see Figure 1).

Figure 1: Plain line: the real path X, dashed lines: virtual paths X(•; µ).

The variation of X is defined as

δX = (δX 1 , δX 2 , . . . , δX N ) := ∂X ∂µ (•; µ) |µ=0 ∈ R 3N .
For each value of µ ∈ (-µ 0 , µ 0 ), we may compute A t0,t1 [X (•; µ)] by replacing X (t) in (2.1) with X (t; µ). The variation of the Hamiltonian action with respect to X between t 0 and t 1 is then defined by

δ X A t0,t1 := d dµ A t0,t1 [X (•; µ)] |µ=0 .
2.1 Principle of stationary action.

The principle of stationary (or least) action states that the real motion X is a critical point of the Hamiltonian action, that is δ X A t0,t1 = 0.

The whole evolution of the system is actually encoded in the above equality. Indeed, let us observe that

δ X A t0,t1 = t1 t0 N i=1 ∂L ∂X i t, X, Ẋ δX i + ∂L ∂ Ẋi t, X, Ẋ dδX i dt dt = t1 t0 N i=1 ∂L ∂X i t, X, Ẋ - ∂ ∂t ∂L ∂ Ẋi t, X, Ẋ δX i (t) dt.
Since X(•, µ) is arbitrary and δX i (t) could be non-zero at any t ∈ (t 0 , t 1 ), the integrand must vanish and we get the so-called Euler-Lagrange equations:

∂ ∂t ∂L ∂ Ẋi t, X, Ẋ - ∂L ∂X i t, X, Ẋ = 0 for all i ∈ 1, N . (2.2)
Time homogeneity and energy conservation. Let us assume that the Lagrangian is independent of time, i.e. L = L X, Ẋ . Taking the scalar product of (2.2) with Ẋi and summing over i ∈ 1, N , we find that

d dt Ẋi ∂L ∂ Ẋi X, Ẋ -L X, Ẋ = 0. (2.3)
The above relation is known as the energy conservation equation, where the total energy of the system would be defined as

E X, Ẋ := Ẋi ∂L ∂ Ẋi X, Ẋ -L X, Ẋ . (2.4)
In general, the Lagrangian is the difference between the kinetic energy T and the potential energy W , that is

L = T -W.
Usually, the kinetic energy is quadratic with respect to Ẋ, i.e.

T = T X, Ẋ = 1 2 ẊT A(X) Ẋ,
where A(X) is a positive definite matrix and W depends only on the generalized coordinates W = W (X). Replacing these expressions into (2.4), one obtains E = T + W which is the total energy of the system.

Case of constrained systems

From the physical point of view, it may be crucial to take into account additional constraints in the system:

• either a geometric constraint F (t, X, Ẋ) = 0 like the search of geodesics on a surface, or the double pendulum problem with two particles connected by rigid rods;

• or an integral constraint (also called isoperimetric constraints) t1 t0 F (t, X, Ẋ)dt = 0 like the famous queen Dido's problem (find the extremal area for a given perimeter) or the catenary (hanging chain) problem.

In that case, we introduce a Lagrange multiplier λ associated to the constraint and replace in the previous calculations L respectively by

• L + λ(t)F ;
• or L + λF with λ independent of time.

For more details on these subjects one may consult classical references on calculus of variations [Sagan, 1992, Gelfand andFomin, 1963].

Newton's equations of movement

As a simple illustration of the previous considerations, let us show how one can recover Newton's second law. Consider a system of N particles animated by the Lagrangian:

L = N i=1 1 2 m i Ẋi 2 -W (X 1 , X 2 , • • • X N ) , (2.5)
where the potential W is a real valued function. Of course this special form of the Lagrangian is obtained as a consequence of physical considerations: existence of an inertial frame, Galileo's relativity principle, time homogeneity etc., see the first chapter of [Landau and Lifschitz, 1976] for more details. In this case we obtain the following Euler-Lagrange equations:

m i Ẍi = - ∂W ∂X i for i ∈ 1, N .
(2.6)

For any i ∈ 1, N , the quantity -∂W ∂X i can interpreted as the force that acts on the i th particle and the Euler-Lagrange equation (2.5) is nothing else but Newton's second law.

Hamilton's principle for continuous deformable media

In this section we propose a formulation of the principle of stationary action in the context of continuous media and derive the corresponding Euler-Lagrange equations. All the notions (virtual displacement, action and variation of the action) introduced in the previous section are generalized to the continuous-media case. Although we use the same notations X and X as in Section 2, we warn the reader that the mathematical objects they describe are of a different nature. This choice is motivated by the fact that in this way, the reader may easily draw parallels and compare the classical particle-setting with the continuous setting.

Continuous deformable media

From a mathematical point of view, a continuous deformable medium (abv. CDM below ) is a pair (Ω 0 , X) such that

• Ω 0 ⊂ R n is a given open set which represents the initial configuration of the medium;

• X : [0, ∞) × Ω 0 → R n is a continuous function satisfying -for all x ∈ Ω 0 , X (0, x) = x; -for all t ≥ 0, X t = X (t, •) : Ω 0 → X t, Ω 0 is a homeomorphism while X t = X (t, •) : Ω 0 → X (t, Ω 0 ) is a diffeomorphism; -for all t ≥ 0, det DX (t, x) > 0.
Physically speaking, we are given a material occupying the initial configuration Ω 0 ⊂ R n , each point x ∈ Ω 0 being occupied by a particle of the medium which will change its position as time evolves. The value X (t, x) represents the position at time t of the particle that was located in x at time 0. The image of Ω 0 through X (t, •), Ω t := X t, Ω 0 = X (t, x) : x ∈ Ω 0 represents the configuration of the medium at time t ≥ 0 (see Figure 2).

Figure 2: Formal representation of the application X t and its image. The domain Ω 0 is seen as the reference space for the motion.

Lagrangian and Eulerian standpoints. There are two different ways of thinking any physical quantity B that characterizes the deformable medium (Ω 0 , X):

• the Lagrangian point of view: B is attached to the particle motion. In this setup B (t, x)

represents the measure at time t of the physical quantity associated to the particle that leaves x at time 0. We will emphasize the Lagrangian character explicitly by an upper script: B . It should be understood automatically that the domain of definition of B is R + × Ω 0 ;

• the Eulerian point of view: B is measured as time evolves, say between t 1 , t 2 in a geometrical point y ∈ t∈[t1,t2] X t, Ω 0 . We will explicitly emphasize the Eulerian character by an e upper script: B e . It should be understood automatically that the domain of definition of

B e is t≥0 {t} × X t, Ω 0 ;
The two points of view are equivalent in the sense that any "Lagrangian physical quantity" can be regarded as an "Eulerian physical quantity" via the family of diffeomorphisms (X (t, •)) t≥0 . More precisely, given B we have

B e (t, y) = B t, X -1 (t, y) for all t ≥ 0 and y ∈ X t, Ω 0 , (3.1)
and conversely

B (t, x) = B e (t, X (t, x)) for all t ≥ 0 and x ∈ Ω 0 . (3.2)
As it will become clearer quite soon, the Lagrangian point of view is appropriate in order to formulate a Hamilton-type principle for CDMs that is close to the classical case briefly presented in Section 2. Thus, we anticipate that we will first obtain equations verified by the Lagrangian physical quantities characterizing the CDM. As we are also interested in getting the equivalent Eulerian equations, we will need some formulae to pass from one formulation to the other. On this point we refer the reader to the introductory section of these notes.

Using the conventions established by relations (3.1) and (3.2) and taking into consideration (1.6) along with (1.7) we have the following relations trace Dψ (DX)

-1 = Dψ ij (DX) -1 ji = (DX) -1 ji ∂ j ψ i = (div ψ e ) , (3.3) 1 det(DX) div (DX) -1 B det DX = (∇B e ) , (3.4)
In particular, taking B = 1 we find the so-called Piola identity

div (DX) -1 det DX = 0. (3.5)
Particle velocity and material derivative In accordance with classical mechanics, we introduce the velocity v such that

Ẋ (t, x) = v (t, x) = v e (t, X (t, x))
which represents the velocity of the particle located at X (t, x). Thus, we have that

Ẋ (t, x) = v e (t, X (t, x)) , X (0, x) = x.
We are now in the position of stating the time derivative law:

∂ t B = (∂ t B e + v e • ∇B e ) . (3.6)
By a material domain we understand the following collection of subsets of R n that is indexed by time: {X (t, ω)} t≥0 . Physically, of course, X (t, ω) represents the region of the space occupied by the particles that were occupying ω at time t = 0.

For a given ω ⊂ Ω 0 it is of physical interest to see how an Eulerian physical quantity evolves in a material domain {X (t, ω)} t≥0 . First, we observe that

Vol (X (t, ω)) = X(t,ω) dy = ω det DX (t, x) dx.
(3.7)

Hence, det DX is related to the local volume change.

More generally, consider now f e : t≥0 {t} × X t, Ω 0 → R and ω ⊂ Ω 0 and using (3.6) let us first write that

X(t,ω) f e (t, y) dy = ω f l (t, x) det DX (t, x) dx from which we infer d dt X(t,ω) f e dy = ω ∂ t f det DX + f ∂ t det DX dx = ω ∂ t f det DX + f det DX (DX) -1 ji ∂ t (DX) ij dx = ω ∂ t f det DX + f det DX (DX) -1 ji Dv ij dx = ω ∂ t f det DX + f (div v e ) det DX dx
where we have used (3.3) to get the last equality. Hence, using (3.6), we obtain

d dt X(t,ω) f e dy = ω (∂ t f e + div (f e v e )) det DX dx = X(t,ω)
{∂ t f e + div (f e v e )} dy.

In particular for f ≡ 1, we infer that

d dt X(t,ω) dy = ω ∂ t det DX = X(t,ω) div v e dy = ω (div v e ) det DXdx.
Since ω is arbitrarily, we obtain that

∂ t det DX = (div v e ) det DX (3.8)
and thus we recover the famous formula for the Wronskian:

det DX = (det DX) |t=0 exp t 0 (div v e ) dt .

Hamilton's principle

The main difference between continuous mechanics and classical finite-degrees-of-freedom mechanics is that the potential energy of a continuous deformable medium is due to the deformation of the medium [Berdichevsky, 2009, Godunov, 1978, Serrin, 1959]. We use this observation in order to postulate a form of the Lagrangian associated to a CDM. Similarly to the discrete case, we define the action of the CDM (Ω 0 , X) between the times t 0 and t 1 as

A t0,t1 [X] := t1 t0 Ω0
Ldx dt, (3.9)

where the Lagrangian L characterizing the CDM (Ω 0 , X) is a given function of the velocity Ẋ, deformation gradient DX, Lagrange coordinates x and time t. For the description of fluids, it is sufficient to take into account only the dependence on d = detDX characterizing the local volume change (see (3.7)) : L = L Ẋ, d, t, x . In the general case, it depends on the whole right Cauchy-Green deformation tensor C = (DX) T DX (in the case of isotropic solids, it depends only on the invariants of C) [Godunov, 1978, Berdichevsky, 2009]. Also, the dependence on the rate of deformation ∂ t d [START_REF] Miles | Weakly dispersive nonlinear gravity waves[END_REF], [Salmon, 1988], [Salmon, 1998], [Gavrilyuk, 1994], [START_REF] Gavrilyuk | A new form of governing equations of fluids arising from Hamilton's principle[END_REF], [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF], [Gavrilyuk, 2011] and its gradients ∇d [Van der Waals, 1979, Truskinovskii, 1982], [START_REF] Gavrilyuk | A model of a plug-chain system near the thermodynamic critical point: connection with the Korteweg theory of capillarity and modulation equations[END_REF], [START_REF] Casal | Connection between the energy equation and the motion equation in Korteweg theory of capillarity[END_REF], [START_REF] Gavrilyuk | Media with equations of state that depend on derivatives[END_REF], [START_REF] Dell'isola | Nucleation of spherical shell-like interfaces by second gradient theory: numerical simulations[END_REF], [Madelung, 1927], [START_REF] Carles | Madelung, Gross-Pitaevskii and Korteweg[END_REF], [Bresch et al., 2019] can be included. In the present notes, we shall restrict ourselves to the most simple situation of classical fluids or a mixture of classical fluids.

Hamilton's principle of stationary action for CDM Suppose that for any time t ∈ [t 0 , t 1 ], Ω t ⊂ R n represents the configuration at time t of the CDM (Ω 0 , X). We consider a family of virtual motions associated to X(t, •), that is a family {X(t, •, µ)} µ∈(-µ0,µ0) verifying the following the properties:

• for all t ∈ [t 0 , t 1 ] and µ ∈ (-µ 0 , µ 0 ), X(t, •, µ) : Ω 0 → R n is a homeomorphism on its image such that X (t, •; µ) |Ω0 is a diffeomorphism on its image;
• for all t ∈ [t 0 , t 1 ] and x ∈ Ω 0 we have that X (t, x; 0) = X (t, x);

• for all x ∈ Ω 0 , X (t 0 , x; µ) = X (t 0 , x) and X (t 1 , x; µ) = X (t 1 , x) for all µ ∈ (-µ 0 , µ 0 );

• for all t ∈ [t 0 , t 1 ], µ ∈ (-µ 0 , µ 0 ) and x ∈ ∂Ω 0 , X (t, x; µ) = X (t, x).

The variation of the family of virtual motions X is defined as ∀t ∈ [t 0 , t 1 ] and ∀x ∈ Ω 0 : δX (t, x) := ∂ ∂µ [X (t, x; µ)] |µ=0 .

(3.10)

In particular, the variation of the Hamiltonian action A between t 0 and t 1 (defined in (3.9)) with respect to X is

δ X A t0,t1 := d dµ A t0,t1 [X] |µ=0 . (3.11)
Then X is a critical point of the Hamilton action (3.9) if, for any X verifying the above properties,

δ X A = 0.
(3.12)

Remark 3.1 The reader will notice the analogy with the classical finite-degrees-of-freedom case presented in Section 2 in the definition of the family of virtual motions: in the classical case, we have fixed the initial point and the endpoint for all curves describing the particles trajectories, while in the case of a continuum we fix the points at the boundary of the whole 4D time-space domain [t 0 , t 1 ] × Ω 0 .

Remark 3.2 Compared to other definitions of Hamilton's action that can be found in the literature (see for instance [Gavrilyuk, 2011]), it is important to note that adopt here a fully Lagrangian standpoint: we integrate the Lagrangian over the domain Ω 0 and consider the Lagrangian variations of the action.

Let us see how we can use the stationary action principle to obtain the corresponding Euler-Lagrange equations. Writing that

d dµ A t0,t1 [X] = t1 t0 Ω0 ∂L ∂ Ẋ ∂ t ∂X ∂µ + ∂L ∂d ∂ ∂µ (det DX) dx dt = t1 t0 Ω0 - ∂ ∂t ∂L ∂ Ẋ ∂X ∂µ + ∂L ∂d det DX trace (DX) -1 D ∂X ∂µ dx dt = - t1 t0 Ω0 ∂ ∂t ∂L ∂ Ẋ + div (DX) -1 ∂L ∂d det DX ∂X ∂µ dx dt (3.13)
and using that X is a critical point of the Hamilton action, we obtain the following Euler-Lagrange equation in Lagrangian coordinates:

∂ ∂t ∂L ∂ Ẋ + div (DX) -1 det DX ∂L ∂d = 0. ( 3
∂ t ∂L ∂ Ẋ Ẋ -L + div (DX) -1 det DX ∂L ∂d Ẋ = 0, (3.15)
Additional state variable. It might also happen that the physical modeling requires the introduction of an extra independent variable (sometimes called order parameter ), say α for the sake of clarity. The state of the system is characterized by the pair (X, α). This should be of course viewed in the Lagrangian setting which can depend on the variable α and its, for example, first time and space derivatives ∂ t α and ∇α (dispersive systems). and adapt Hamilton's principle by imposing that (X, α) is critical. We skip the technical details and just note that the Euler-Lagrange equations in the Lagrangian coordinates in this case become

         ∂L ∂α - ∂ ∂t ∂L ∂(∂ t α) -div ∂L ∂(∇α) = 0, (3.16a) ∂ t ∂L ∂ Ẋ + div (DX) -1 det DX ∂L ∂d = 0. (3.16b)
As before, the energy conservation equation can be derived as a consequence of the governing equations:

∂ t ∂L ∂ Ẋ Ẋ + ∂L ∂(∂ t α) ∂ t α -L + div (DX) -1 det DX ∂L ∂d Ẋ + ∂L ∂(∇α) ∂ t α = 0.
(3.17)

Such an approach was successfully applied for the Eulerian description of wave propagation in fluids containing gas bubbles [Drui et al., 2019, Gavrilyuk andSaurel, 2002], in dispersive hydrodynamics [START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF] and even in quantum mechanics [START_REF] Dhaouadi | Extended Lagrangian approach for the defocusing nonlinear Schrödinger equation[END_REF].

For the sake of clarity, in the following we consider the most simple situation where the dynamics does not depend on the space and time derivatives of α and Eq. (3.16a) reduces to ∂L ∂α = 0.

(3.18)

Mixtures as CDMs with different velocities

The multi-phase flow modeling is a vast and challenging topic. Multiple important applications necessitate the construction of physically reliable and mathematically well-posed models. Classical averaging approach [START_REF] Drew | Theory of multicomponent fluids[END_REF], [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF], [Nigmatulin, 1990] usually produces an underdetermined system of equations. Different closure hypotheses, even if they are reasonable, can drastically change the mathematical structure of the governing equations.

A simple idea could be to formulate the Lagrangian of a complex multi-phase system directly in terms of average variables. At least it allows us to keep the same mathematical structure of the governing equations. In a second step, one takes in consideration friction forces which, at a mathematical level, amounts to add dissipative terms that are compatible with the entropy inequality. Thus, the main issue is the formulation of the non-dissipative equations (Euler-Lagrange equations) defined solely by the Lagrangian of the system. We want to model the case of a continuous deformable medium that is composed of two types of particles. In the most general case, each component is animated by its own movement. We consider the two CDM (Ω 0 , X 1 ), (Ω 0 , X 2 ) with the property that

X 1 t, Ω 0 = X 2 t, Ω 0 ∀ t ≥ 0. (3.19)
We are saying that at t = 0, any position x ∈ Ω 0 is occupied by two particles of different types. As time evolves, the particles change positions but in such a way that at any t > 0, the global domain occupied by the first specie is the same as the domain occupied by the second specie. Hence, no pure phase is allowed and (3.19) can be seen as a "no-segregation property". Note also that at time t a geometrical point y ∈ X 1 t, Ω 0 = X 2 t, Ω 0 is not necessarily occupied by particles that left the same position at the initial time (see Figure 3).

Figure 3: Two formal representations of the system. On the left, the initial configuration of the system Ω 0 and its image Ω t = X 1 (t, Ω 0 ) = X 2 (t, Ω 0 ). On the right, the two reference spaces

Ω 1 = X -1 1 (t, Ω t ), Ω 2 = X -1 2 (t, Ω t )
. One passes from one reference space Ω j to the other one Ω i through the application X

-1 i • X j (t, •).
What remains for a complete description of the mixture is to say "how much volume at a point" does a particle of type a ∈ {1, 2} occupy. To that end, we introduce α a (t, x) = the volume fraction at time t of the particle of type a measured at X a (t, x) .

More precisely, for a given point (t, x), we denote y = X a (t, x) and consider a representative mixture volume containing y, ω = ω 1 ∪ ω 2 where ω a is the volume occupied by the phase a. We define the characteristic function of phase a as

χ e a (t, z) = 1 if z ∈ ω a , 0 otherwise.
Then α e a (t, y) := ω χ e a (t, z) dz

Vol(ω) = Vol(ω a ) Vol(ω) .
In particular, with hypothesis (3.19) we impose 0 < α e a (t, y) < 1 ∀ t ≥ 0, y ∈ X a (t, Ω 0 ).

Remark 3.3 It is important to note that the definition of the volume fraction depends on the averaging process which is performed. Several averaging methods (e.g. ensemble averaging, statistical averaging, time averaging, etc.) may be found in the literature and we refer the interested reader to the books [START_REF] Drew | Theory of multicomponent fluids[END_REF], [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF]] and [Nigmatulin, 1990].

The two volume fractions are linked through the relation

α 1 (t, x) + α 2 t, X -1 2 (t, X 1 (t, x)) = α 1 t, X -1 1 (t, X 2 (t, x)) + α 2 (t, x) = 1. (3.20)
Indeed, the sum of the volume fractions is 1 at the geometrical point y = X 1 (t, x), position which is occupied at time t by the particle of the second specie that left X -1 2 • X 1 (t, x) at time t 0 = 0. In order to ease the reading, when there is no possible confusion, we shall drop the upper-script in the computations that follow.

We would like to emphasize the fact that the choice of a physically meaningful Lagrangian L mix is always a difficult problem. A priori, it cannot be reduced to a linear combination of the pure phases Lagrangians because of the interaction between phases. In particular, a new scalar variable | Ẋ1 -Ẋ2 | should be added for the description of the virtual mass effects (when a moving particle of type 1 will displace surrounding fluid particles of type 2). The Lagrangian can also depend on topological characteristics of multiphase flows, as, for example, on the total area of contact between two fluids, i.e. the interfacial area (see the recent study [START_REF] Cordesse | Derivation of a two-phase flow model with two-scale kinematics and surface tension by means of variational calculus[END_REF]). In the present notes, we shall stick to most simple physics and assume that the Lagrangian characterizing the dynamics of the mixture has the form

L mix = L 1 Ẋ1 , detDX 1 , α 1 , x + L 2 Ẋ2 , detDX 2 , α 2 , x .
The associated Hamiltonian action between the times t 0 and t 1 is then

A t0,t1 [X 1 , X 2 , α 1 ] = A 1;t0,t1 [X 1 , α 1 ] + A 2;t0,t1 [X 2 , α 2 ] (3.21) = t1 t0 Ω0 L 1 Ẋ1 , det DX 1 , α 1 , x + L 2 Ẋ2 , det DX 2 , α 2 , x dx dt.
We postulate that the equations governing (X 1 , X 2 , α 1 ) are obtained as the Euler-Lagrange equations associated to the action A. These equations are obtained by imposing that the variations of A with respect to X 1 , X 2 and α 1 are all 0. In the derivation of the equations it is important, even crucial, to remember the relation (3.20) which implies a cross-contribution of L 2 in δX 1 , and conversely.

We begin with the variation of A with respect to α 1 . Of course, the delicate part is to obtain the contribution coming from A 2 . Owing to the relation (3.20), we put A 2 under the form

A 2;t0,t1 = t1 t0 Ω0 L 2 Ẋ2 , det DX 2 , α 2 , x dx dt = t1 t0 Ω0 L 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1(t,x ) det D(X -1 2 • X 1 ) |(t,x ) dx dt = t1 t0 Ω0 L 2 Ẋ2|X -1 2 •X1(t,x ) , det DX 2|X -1 2 •X1(t,x ) , 1 -α 1 (t, x ), (X -1 2 • X 1 )(t, x ) × det D(X -1 2 • X 1 ) |(t,x ) dx dt.
Thus we get that

δ α1 A t0,t1 [X 1 , X 2 , α 1 ] = t1 t0 Ω0 ∂L 1 ∂α 1 Ẋ1 , det DX 1 , α 1 , x - ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 det D X -1 2 • X 1 dx dt,
and thus we obtain the first local equation

∂L 1 ∂α 1 Ẋ1 , det DX 1 , α 1 , x = ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 det D X -1 2 • X 1 . (3.22)
Next, we study the variation with respect to X 1 . Again, the variation of A 1;t0,t1 is treated as in (3.13), we focus only on the contribution coming from A 2;t0,t1 :

δ X1 A 2 = t1 t0 Ω0 δ X1 L 2 Ẋ2 , det DX 2 , α 2 , x dx dt = t1 t0 Ω0 ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x δ X1 α 2 dx dt. But writing that α 2 (t, x) = 1 -α 1 t, X -1 1 (t, X 2 (t, x)) , we see that δ X1 α 2 (t, x) = -Dα 1|(t,X -1 1 •X2(x)) δ X1 (X 1 ) -1 (t, X 2 (t, x)) . We get then δ X1 α 2|(t,X -1 2 •X1(t,x)) = -Dα 1 δ X1 (X 1 ) -1 (t, X 1 (t, x)) = Dα 1 D(X -1 1 ) |(t,X1(t,x)) δX 1 (t, x).
The passage to the last equation is obtained by differentiating the relation X

-1 1 • X 1 = Id with respect to µ, i.e. ∂X -1 1 ∂µ (t, X 1 (t, x)) + D(X -1 1 ) |(t,X1) ∂X 1 ∂µ (t, x) = 0,
and for µ = 0 we get

δ X1 (X 1 ) -1 (t, X 1 (t, x)) = -D(X -1 1 ) |(t,X1) δX 1 (t, x) .
We may then write that

δ X1 A 2;t0,t1 = t1 t0 Ω0 ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x δ X1 α 2 dx dt = t1 t0 Ω0 ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x δ X1 α 2 |X -1 2 •X1 det D X -1 2 • X 1 dx dt = t1 t0 Ω0 ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 Dα 1 (DX 1 ) -1 det D X -1 2 • X 1 δX 1 dx dt.
Consequently, we have

δ X1 A t0,t1 = - t1 t0 Ω0 ∂ ∂t ∂L 1 ∂ Ẋ1 + div (DX 1 ) -1 ∂L 1 ∂d 1 det DX 1 δX 1 dx dt + t1 t0 Ω0 ∂L 2 ∂α 2 |X -1 2 •X1 Dα 1 (DX 1 ) -1 det D X -1 2 • X 1 δX 1 dx dt,
which yields the Euler-Lagrange equation associated to the variation in X 1 :

∂ ∂t ∂L 1 ∂ Ẋ1 + div (DX 1 ) -1 det DX 1 ∂L 1 ∂d 1 = ∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 Dα 1 (DX 1 ) -1 det D X -1 2 • X 1 . (3.23)
Obviously we obtain a similar equation when considering the variation with respect to X 2 :

∂ ∂t ∂L 2 ∂ Ẋ2 + div (DX 2 ) -1 det DX 2 ∂L 2 ∂d 2 = ∂L 1 ∂α 1 Ẋ1 , det DX 1 , α 1 , x |X -1 1 •X2 Dα 2 (DX 2 ) -1 det D X -1 1 • X 2 . (3.24)
Thus the Euler-Lagrange equations associated to the variation of the Hamiltonian action are given by (3.22), (3.23) and (3.24). As it will be seen in application, the Eulerian version of the above equations have a simpler form.

Recovering the Euler equations

Formally, a fluid is a CDM (Ω 0 , X) to which we associate two state functions ρ , η representing respectively the density and the entropy of the fluid. We denote

ρ |t=0 = ρ 0 , η |t=0 = η 0
the density and the entropy at the initial time t 0 = 0. In the calculations that follow it may be sometimes convenient to work with the inverse of the density,

τ := 1 ρ , (4.1)
which is called the specific volume.

Before derivating the equations describing the fluid motion, we need to precise two fundamental assumptions.

Mass conservation

We assume that the physical system is closed, i.e. no matter is created or destroyed in the system. Then the mass M e (t, X(t, ω)) at time t ≥ 0 of the material volume X (t, ω) defined as

M e (t, X(t, ω)) := X(t,ω) ρ e (t, y) dy = ω ρ (t, x) det DX (t, x) dx, remains constant in time. That is, for all ω ⊂ Ω 0 and t ≥ 0 ω ρ (t, x) det DX (t, x) dx = ω ρ (0, x) det DX (0, x) dx = ω ρ 0 (x) dx, which amounts to ask that ρ (t, x) det DX (t, x) = ρ 0 (x) ∀ t ≥ 0, x ∈ Ω 0 . (4.2)
Obviously the last identity, along with (4.1), yields

τ (t, x) = det DX (t, x) ρ 0 (x) .
Transport of the entropy We assume that the flow is reversible, i.e. the entropy is conserved along particle trajectories

η (t, x) = η 0 (x) ∀ t ≥ 0, x ∈ Ω 0 . (4.3)

The case of a compressible fluid

We define the Lagrangian kinetic energy density of the fluid as

T Ẋ, x := 1 2 ρ 0 (x) | Ẋ| 2 , (4.4)
and the Lagrangian potential energy density

W (det DX, x) := ρ 0 (x) ε det DX ρ 0 (x) , η 0 (x) = ρ 0 (x) ε τ , η . (4.5)
where we call ε = ε(τ, η) the (Lagrangian) internal energy per unit mass which related to the (Lagrangian) local pressure through the relation:

p (τ, η) = -∂ τ ε (τ, η) . (4.6)
We postulate that the Lagrangian that encodes the equations of movement of the fluid is defined as:

L Ẋ, det DX, x := T Ẋ, x -W (det DX, x) = 1 2 ρ 0 (x) | Ẋ| 2 -ρ 0 (x) ε det DX ρ 0 (x) , η 0 (x) .
We immediately see that ∂L ∂ Ẋ Ẋ, det DX, x = ρ 0 (x) Ẋ, while, using (4.2), we get that

∂L ∂d Ẋ, det DX, x = -ρ 0 ∂ τ ε det DX ρ 0 , η 0 1 ρ 0 = p τ , η .
Then, recalling that Ẋ(t, x) = v (t, x), the Euler-Lagrange equation (3.14) becomes

ρ 0 ∂ t v + div (DX) -1 p τ , η det DX = 0.
Multiplying the last equation by 1/ det DX and using the lagrangian-eulerian passing formulae (3.2), (3.6) and (3.4) we get that ρ e (∂ t v e + v e • ∇v e ) + ∇p (τ e , η e ) = 0.

Let us summarize the above discussion. The Lagrangian equations for the evolution of a fluid are given by

   ρ det DX = ρ 0 , η = η 0 , ρ 0 ∂ t v + div (DX) -1 p τ , η det DX = 0.
In Eulerian coordinates, we recover the compressible Euler system

   ∂ t ρ e + div (ρ e v e ) = 0, ∂ t η e + v e
• ∇η e = 0, ρ e (∂ t v e + v e • ∇v e ) + ∇p (τ e , η e ) = 0.

(4.7)

Incompressible Euler equations

When the volume of any material element (X (t, ω)) t≥0 does not change in time, the flow is said to be incompressible. Remembering formula (3.7), the incompressibility condition reads

det DX (t, x) = 1 ∀ t ≥ 0, x ∈ Ω 0 , (4.8)
from which we infer

ρ (t, x) = ρ 0 (x) ∀ t ≥ 0, x ∈ Ω 0 .
Since det DX is constant, there is no energy coming from the deformation of the fluid and a naive choice for the Lagrangian would be:

L Ẋ, det DX, x = T ( Ẋ, x) = 1 2 ρ 0 (x)| Ẋ| 2 .
Nonetheless, recalling the discussion of Subsection 2.2, the constraint (4.8) leads us to consider the modified Lagrangian L Ẋ, det DX, p := 1 2 ρ 0 (x) | Ẋ| 2 + p (det DX -1) .

(4.9)

The new variable p is called pressure of the fluid. However, it is important to note that the two pressures p and the "compressible pressure" (4.6) have different meanings. For an incompressible flow, the pressure is seen as a Lagrange multiplier associated to the constraint det DX (t, x) = 1.

Considering the variations with respect to Ẋ and p of the modified Lagrangian, we have

       ∂ L ∂ Ẋ = ρ 0 (x) Ẋ, ∂ L ∂d = p ,
so that the Euler-Lagrange equations now read

       ρ = ρ 0 , η = η 0 , ρ 0 ∂ t v + div (DX) -1 p = 0, det DX = 1.
Using (3.2), (3.6) and (3.4), we easily translate the first three equations in the Eulerian variables.

In the Eulerian setting, we recover from the last equation the famous relation div v e = 0 (we recall (3.8)). We thus have

       ∂ t ρ e + u e
• ∇ρ e = 0, ∂ t η e + u e • ∇η e = 0, ρ e (∂ t v e + v e • ∇v e ) + ∇p e = 0, div v e = 0.

Dynamics of two-phase flows

We have previously explained in Section 3.3 what we mean by a mixture of two species a ∈ {1, 2} occupying an initial configuration Ω 0 and that are animated by the two motions X 1 , X 2 . Under the "no-segregation assumption" (3.19), we saw that the system is also described by volume fractions of the two components which, as before, we call α 1 , α 2 . The relation linking them is

α 1 (t, x) + α 2 t, X -1 2 • X 1 (t, x) = 1, (5.1)
where, by X -1 2 •X 1 (t, x), we should understand X -1 2 (t, X 1 (t, x)). Each component of the mixture is characterized by the state variables ρ a and η a . We define the apparent density ρa of the acomponent which measures the density of the a-component in the mixture. It is different from the characteristic density, both being linked through the relation

ρ a = ρ a α a for a ∈ {1, 2} .
We emphasize that the quantities ρ a , η a follow the particles at their own characteristic velocity so that ρ 1 (t, x) , ρ 2 (t, x) do not represent in general the density in a same physical point ! As in the single-phase case presented in the previous section, we assume that the mass of each phase is conserved, that is ρ a det DX a = ρa,0 for a ∈ {1, 2} , (

where ρ1,0 , ρ2,0 : Ω 0 → [0, ∞) are the apparent densities at the initial time. Moreover, we assume that the entropies of the constituents remain constant along particle paths meaning that η a = η a,0 , (5.3)

where η 1,0 , η 2,0 : Ω 0 → R are the entropies at the initial time. We consider two functions modeling the internal energies of the components

ε a = ε a τ a , η a ,
where τ a = 1/ρ a and we denote by p a τ a , η a = -∂ τ a ε a τ a , η a the corresponding pressures.

One velocity two-phase models

A first class of models for mixtures is obtained by assuming that the particles move with the same velocity, i.e. v 1 = v 2 = v, which amounts to making the following simplification:

X 1 = X 2 = X.
Therefore the following equations hold true

α 1 (t, x) + α 2 (t, x) = 1 ∀ t ≥ 0, x ∈ Ω 0 ,
ρ a det DX = ρa,0 , η a = η a,0 for a ∈ {1, 2}.

We define the Lagrangian kinetic energy density (resp. Lagrangian potential energy density) of the mixture as the sum of the kinetic (resp. potential) energies of each component (of course this choice is questionable, see the discussion in Subsection 3.3):

T Ẋ, x = 1 2 (ρ 1,0 (x) + ρ2,0 (x)) | Ẋ| 2 , W (det DX, α 1 , x) = ρ1,0 (x) ε 1 α 1 det DX ρ1,0 (x)
, η 1,0 (x) + ρ2,0 (x) ε 2 α 2 det DX ρ2,0 (x) , η 2,0 (x) .

We postulate that the Lagrangian which encodes the dynamics of the mixture is defined as

L Ẋ, det DX, α 1 , x = T Ẋ, x -W (det DX, α 1 , x).
The Euler-Lagrange equation associated to the variation of X is computed following the same steps as in the previous sections:

∂L ∂ Ẋ Ẋ, det DX, α 1 , • = (ρ 1,0 + ρ2,0 ) Ẋ = (ρ 1,0 + ρ2,0 ) v , and ∂L ∂d Ẋ, det DX, α 1 , • = α 1 p 1 + α 2 p 2 =: p ,
where p is defined as the pressure of the mixture. Consequently, the Euler-Lagrange equation reads (ρ 1,0 + ρ2,0 ) ∂ t v + div (DX) -1 p det DX = 0.

The system is closed by considering the variation according to the volume fraction which leads to (3.18) (cf Section 3.2), i.e. ∂L ∂α 1 ( Ẋ, det DX, α 1 , x) = 0.

(5.4)

Since the kinetic energy is independent of α 1 , we just have to analyze the variation of the potential energy with respect to α 1 . We observe that ∂W a ∂α a = ρa,0 ∂ τ ε a α a det DX ρa,0 , η a,0 det DX ρa,0 for a ∈ {1, 2}.

Hence, taking into account the constraint

α 1 + α 2 = 1, Equation (5.4) yields p 1 det DX = p 2 det DX.
Gathering the equations, we obtain that

         η a = η a,0 , α a ρ a det DX = ρa,0 , ρ 0 ∂ t v + div (DX) -1 p det DX = 0, p = p 1 = p 2 ,
where have set ρ 0 := ρ1,0 + ρ2,0 .

Using the relations (3.4), (3.3) and (3.6), we obtain the corresponding Eulerian equations which are

       ∂ t η e a + v e • ∇η e a = 0, ∂ t (α e
a ρ e a ) + div (α e a ρ e a v e ) = 0, ρ e (∂ t v e + v e • ∇v e ) + ∇p e = 0, p e = p e 1 = p e 2 ,

(5.5)

where the mixture density is defined as

ρ e := α e 1 ρ e 1 + α e 2 ρ e 2 .
The system (5.5) can obviously be generalized to any number of components. It is sometime called Kapila's model and was obtained in [START_REF] Kapila | Twophase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF].

From the mathematical point of view, one can show that the equations (5.5) are hyperbolic if ∂p a ∂τ a < 0 for any a ∈ {1, 2}.

(5.6)

The sound speed c W in such a system is given by the Wood formula :

1 ρ e c 2 W = a α a ρ e a c 2 a , (5.7)
where c a is the sound velocity of a -th component, i.e.

c 2 a = -τ 2 a ∂p a ∂τ a .

Two velocities, one pressure models

We assume that the local kinetic and potential energy densities are given by

     E a Ẋa , x = 1 2 ρa,0 (x) | Ẋa | 2 , W a Ẋa , det DX a , α a , x = ρa,0 ε a α a det DX a ρa,0 , η a,0 .
As we discussed in Section 3.3, the Lagrangian encoding the dynamic of the mixture is given by

L = L 1 + L 2 = (T 1 -W 1 ) + (T 2 -W 2 )
where, of course, it is to be considered that X 1 , X 2 and α 1 are independent. In order to obtain the Euler-Lagrange equations we just have to express the abstract equations (3.22), (3.23) and (3.24) obtained in Section 3.3. We begin with ∂L a ∂α a = -∂W a ∂α a = ρa,0 p a α a det DX a ρa,0 , η a,0 det DX a ρa,0 = p a τ a , η a det DX a .

Also, we have that

∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 det D X -1 2 • X 1 = p a 1 -α 1 ρ 2 • X -1 2 • X 1 , η a det DX 2 X -1 2 • X 1 det D X -1 2 • X 1 = p a 1 -α 1 ρ 2 • X -1 2 • X 1 , η a • X -1 2 • X 1 det DX 1 .
(5.8)

We see that the equation becomes

p 1 α 1 (t, x) ρ 1 (t, x) , η 1 (t, x) det DX 1 = p 2 1 -α 1 ρ 2 • X -1 2 • X 1 , η 2 t, X -1 2 • X 1 (t, x) det DX 1 (5.9)
Obviously the last equation says that the Eulerian pressures are equal p 1 (τ e 1 , η e 1 ) = p 2 (τ e 2 , η e 2 ) .

(5.10)

Let us express the equation coming from the variation of X 1 . Using (5.8) along with (5.9) we get that

∂L 2 ∂α 2 Ẋ2 , det DX 2 , α 2 , x |X -1 2 •X1 Dα 1 (DX 1 ) -1 det D X -1 2 • X 1 = p 2 1 -α 1 ρ 2 • X -1 2 • X 1 , η 2 • X -1 2 • X 1 det DX 1 Dα 1 (DX 1 ) -1 = p 1 det DX 1 Dα 1 (DX 1 ) -1
and thus the corresponding Euler-Lagrange equation is Gathering all the equations, we find that (5.12)

ρa,0 ∂ t v a + div (DX a ) -1 α a p a det DX a = det DX a p 1 Dα a (DX a ) -1 (
               η a = η a,0 , α a ρ a det DX a = ρα,0 , ρa,0 ∂ t v a + div (DX a ) -1 α a p α det DX a = det DX a p a Dα a (DX a ) -1 , α 1 + α 2|X -1 2 •X1 = 1, p 1 = [p 2 ] |X -1 2 •X1
A simple analysis shows that, even if the equations admit the conservation of the total momentum and total energy, they cannot be rewritten as a conservative system of equations: the number of conservation laws admitted by the system is smaller than the number of unknowns. A non-classical definition of weak solutions is thus needed. Moreover, for small relative velocity w = v 1 -v 2 the equations are not hyperbolic: complex characteristics appear. Thus the well-posedness of the Cauchy problem for (5.12) is questionable.

A possibility to evacuate the problem of complex characteristics is to use relaxation methods. Let us consider the system (5.13)

                                  
Here θ a are the phase temperatures : θ a = ∂ ηa ε a , p I and v I are the interface pressure and velocity, respectively, λ is the friction coefficient, and f a are the entropy production terms. The parameter µ is large, so we formally recover the "equal pressure" model as µ → ∞ whereas we recover a "one-velocity" model in the limit λ → ∞.

The total energy equation and the following entropy inequality allow us to find simple (phenomenological) expressions for f a [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF], Gavrilyuk, 2011].

The system (5.13) is also called Baer-Nunziato model [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] (in the case v I = v 1 , p I = p 2 or vice versa). It is weakly hyperbolic, i.e. all eigenvalues are real but one cannot guarantee a full set of the eigenvectors in the resonant case (eigenvectors can coincide on some hypersurfaces [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF], [START_REF] Andrianov | The Riemann problem for the Baer-Nunziato two-phase flow model[END_REF]). The choice of interfacial variables v I and p I is not unique see [START_REF] Hérard | A simple method to compute standard two-fluid models[END_REF] and [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF].

Let us mention that the relaxation equation for the volume fraction (last equation of (5.13)) can be obtained from the barotropic Navier-Stokes equations by a homogenization procedure [START_REF] Bresch | A multi-fluid compressible system as the limit of weak solutions of the isentropic compressible Navier-Stokes equations[END_REF], Bresch and Hillairet, 2019, Bresch and Hillairet, 2015, Hillairet, 2007, Bresch et al., 2020]. Loosely speaking, if the density of a fluid mixture wildly oscillates between two reference densities, the volume fractions α e a can be interpreted as the probability to find a particle of fluid a at a given position and time instant. The coefficient µ in such a homogenization procedure is the inverse of the mean viscosity of the two phases.

Conclusion

The aim of this short course was the introduction into variational approach to the modeling of fluid mixtures. The only scalar function we needed to know is the full energy of the system (or, equivalently, the Lagrangian of the mechanical system). Due to the Noether theorem, this approach guarantees the fulfillment of the basic physical conservation laws, but, a priori, not wellposedness in the sense of Hadamard of the governing equations. The introduction of dissipative terms can provide the well-posedness. The simplest and most popular models of fluid mixtures (5.5) and (5.13) were derived.

The models obtained can also be used for the description of the interface between pure fluids. For example, in the problem of shock-interface interaction the interface between fluids becomes unstable (Meshkov-Richtmyer instability) and can be considered as a mixture of fluids [Karni, 1994, Saurel and[START_REF] Saurel | [END_REF]. This method called also the method of diffuse interfaces, was recently extended to the problems involving solid-fluid mixtures [START_REF] Ndanou | Multi-solid and multifluid diffuse interface model: Applications to dynamic fracture and fragmentation[END_REF] and to the description of interfaces with capillary effects [START_REF] Schmidmayer | A model and numerical method for compressible flows with capillary effects[END_REF].

The method can be extended to more complex situations involving mixtures of continua depending on higher order deformation gradients and the time derivatives of the deformation gradient, like bubbly fluids or the Green-Naghdi equations for dispersive shallow water flows (see [Gavrilyuk, 2011]).

  .14) Similarly to the discrete case (cf (2.3)), when the Lagrangian does not depend explicitly on time, L = L Ẋ, d, x , one derives an energy conservation equation by combining (3.3), (3.8) and (3.14):

  2 (∂ t v e 2 + v e 2 • ∇v e 2 ) + ∇(α e 2 p e 2 ) = α e 2 ∇p e I -λ(v e 2 -v e 1 ),α e 1 + α e 2 = 1, ∂ t α e 1 + v e I • ∇α e 1 = µ(p e 1 -p e 2 ).

  .

	Using the relations (3.4) , (3.3) and (3.6) we obtain the corresponding Eulerian equations which
	are	          	∂ t η e a + v e a • ∇η e a = 0, ∂ t (α e a ρ e a ) + div (α e a ρ e a v e a ) = 0, ∂ t (α e a ρ e a v e a ) + div (α e a ρ e a v e a ⊗ v e a ) + ∇(α e a p e a ) = α e a ∇p e a , α e 1 + α e 2 = 1, p e 1 = p e 2 .


Luc Leviandier

W Alpers

E Barthelemy

P Brandt

P V Brennan

Glodina Connan

A C Edwards

René Garello

Josep A Rodenas Navas

H Griffiths

J Rodenas

D Renouard

A Rubino

J O Thomas

K Woodbridge

MORSE: mesoscale ocean radar signature experiments

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

It is now well established that satellite-borne SAR (synthetic aperture radar) is a powerfull mean to image those ocean surface phenomena like swells or surface slicks that modify its roughness. Besides, internal dynamics of the ocean (internal waves and tides, mesoscale oceanographic features) interacts with its surface roughness through associated surface currents. Among these, long internal waves generated by tidal flow over the slope and propagating on the continental shelf are very energetic and have a considerable oceanographic and economic impact in the coastal zones. So local or global wave information is needed by a wide range of users in various activities like off-shore drilling and mining, ship routing, fisheries and aquaculture, environmental monitoring, acoustic communication, acoustic detection or coastal engineering in addition to physical oceanography itself. Nevertheless, at the present time, radar imagery is not widely used by the potential users community. This is due to the lack • of sufficiently reliable models able to retrieve quantitative results on the 3D structure of the hydrodynamic processes,

• of user-friendly software tools able to automatically process the images, releasing the user from a complex "manual" interpretation task.

The objective of the project is then to furthering the understanding of the physical processes involved in radar signatures of internal waves in the ocean. Tools are developped that will contribute to fill the gap between availability of images and routine use of them.

MORSE is a multidisciplinary approach in the above mentionned direction. It combines:

• Laboratory experiments using a dedicated small W-band (3 mm wavelength) radar imaging surface wave fields produced in laboratory tank and flume in Grenoble.

•

Theoretical research aimed at assessing and developping a three scale hydrodynamic model to improve radar backscattering predictions. The project is structured into four inter-linked technical and scientific parcels, namely:

• the development of image processing algorithms for automatic extraction of surface wave parameters from radar images (together with their implementation and testing on ERS 7 -l and 2 imagery).

•

Experimental and theoretical research on basic non linear hydrodynamic interactions between internal waves and surface gravity waves/capillary waves.

•

The laboratory W-band SAR design, construction and test for tank experiments.

• Research on radar backscatter modeling (including modelling, especially for C-band, and inversion).

The following sections of this paper outline characteristic results obtained in these areas.

With regards to software products, three of the above listed areas, the two first and the last one will result in The MORSE consortium comprises four university group (UCL, IfM, ENSTB, LEGI) and two commercial companies (OCS, TMS)

1997)).

The second part of this section is devoted to another approach for detection and location of internal waves using wavelets.

A) WASP

WASP is the practical wing of the overall software approach and is based on past OCS experience and ENS TB' s state-of the-art knowledge of wave parameter extraction. WASP incorporates an image handling segment (images readable from CD-Rom, tape or hard disk) which includes a bespoke (and simple) macro-language to enable users to define any new image data structures (within reason) so that such images may be read in for further processing. This should future-proof the system to some extent Sampling, clipping, selection, and sub-image export and storage operations are allowed.

An image catalogue/index system is being included in order to supply a complete image/sub image/process audit. In the next segment, various high spatial-resolution parameter extraction algorithms (supplied by ENSTB) may be applied to sub-images to obtain a map of dominant and sub-dominant wavevectors over the image.

Numerical data may be exported. The numerical data may be used to generate an image-overlay curvi-linear coordinate system.

A third component of the system will enable 'transects' (grey level profiles along specified paths) to be obtained and stored.

B) Detection and Location of Internal Wave from SAR Images using wavelets

ENSTB is currently developing a prototype system for automatic detection, verification, assesment and characterisation of internal waves from satellite images. The long-range objective of this work is to develop image and signal processing algorithms that would allow a computer to produce analysis product as given by human analyst given satellite imagery. This is a difficult image analysis and understanding problem due to a lack of precise mathematical descriptions of the ocean features, coupled with the time-varying nature of these features and the compilation of frequent partial obscuration by cloud cover. The total automation of this task is not presently feasible, but efforts toward that end will result in many tools that the analyst can use to reduce labor requirements and improve the quality and consistency of the analysis. Results so far indicate that it should be possible to reach this goal (R6denas et al., 1997).

Here, we present a case study for the detection and location of internal wave patterns from high-resolution SAR images collected either by the European satellites ERS-1/2 or the Canadian RADARSA T over the Strait of Gibraltar or the Andaman Sea.

I) Automatic Internal Wave Detection and Location

In Fig. I, the main modules of the internal wave detector and location system are shown. Our strategy is to detect all possible internal waves, and then locate each individual internal wave packet. Methods for mesoscale feature extraction from SAR images often apply some kind of filtering as a first step. Therefore, the first step in our system involves averaging the original SAR image, 8000 ;x 8000 pixels in size, using a 8 x 8 pixel nonoverlapping window. This produces an image 1/8 the size of the original, thus reducing subsequent processing times. The image also exhibits substancially reduced speckle noise.

In this non-supervised internal wave detection approach, we begin by separating land from sea in order to avoid false alarms in the detection process. Extraction of the shoreline in SAR images is a difficult task to perform using simple image processing operations such as grey-value thresholding, due to the presence of speckle and because the signal returned from the sea surface may be similar to that from the land. Thus, an automatic method for detecting the shoreline accurately and efficiently from SAR images is also introduced in the system. It would be possible to replace the coastline detection step by a manual delineation of the coastline obtained by digitising a map. However, changes that have evolved since the map was produced, such as small sandbanks near the coastline, may then be missed. Then, the image is wavelet transformed with various scales, based on the Wavelet Transform Modulus Maxima (WTMM) method. The wavelet transform gives spectral decomposition via the scale concept. By decomposing signals into elementary building blocks that are well localised both in space and scale, the wavelet transform can characterise the local regularity of signals, such as internal waves, via the study of the wavelet local modulus maxima across decomposition scales. Internal waves features are detected at multiple scales and characterised by their positions and magnitude. Then, a method for internal wave feature enhancement, based on an• edge linking algorithm, is required in order to link these features. At this stage, we also discard significant local maxima belonging to the land class by means of the sea-land mask obtained from the coastline detection method. The remaining features are significant features of the important structures in the SAR image, e.g. internal waves, oil slicks, eddies, natural films, etc. Once we have detected the significant features in the SAR image, some a priori oceanographic knowledge usually is required on the geometric shape and size of the internal wave signatures and look-alikes for the identification and extraction of internal wave features. For this purpose, we use the phase information given by the wavelet maxima for all the significant detected features from the SAR image. Finally, a clustering algorithm, based on the neighboring features, is applied for internal wave packet location in the SAR image.

2) Results and Discussion

In this sub-section, we present several results of applying the proposed method for internal wave detection and location to two SAR images containing inernal wave patterns. The first SAR image to be analysed is show in Fig. 2(a). This is a 3-look ERS-1 SAR image from the Strait of Gibraltar (Spain), 8000 x 8000 pixels in size, with a pixel spacing of 12. 5 m x 12. 5 m, covering an area of approxim a tely 100 km x 100 km and acquired on September 10, 1993, is used as first test image. The image has been displayed using histogram equalisation, to reduce loss of contrast in the reproduction process. This image shows roughness patterns associated with a short wavelength soliton-like internal wave train propagating eastwards into the Strait. It also shows other oceanic processes as for instance : possible oil slicks, ship and ship wakes, natural films, and dark spots due to wind conditions. In the image, these several possible internal wave look-alikes are also underlined. The proposed method for internal wave detection and location has then been tested on this SAR image. From Fig. 2(b ), it can be seen that internal wave features are well detected and located against background noise or other internal wave look-alike phenomena. Mainly six internal wave packets are detected, and located in the SAR image. From the observation of the wavelengths and crests of the detected internal wave packets, we can estimate that the largest packets are propagating westwards, since the spacing between waves increases away from the leading edge, as to be expected in dispersive waves. The image also shows the detected coastline from sea-land classification. Fig. 3(a) shows a RADARSAT SAR image from the Andaman Sea which was acquired on May 4, 1996, covering an area of 300 km x 900 km. This image shows the interactions between several internal soliton packets propagating westwards in the Andaman Sea. The attractive of this image remains in its difficulty to well detect internal wave features due to the high amount of oceanic information in the sea surface. In the image, we can observe Top views of the measurement field at successive times allow to compute the refraction angle between the incident and the transmitted wave over the shelf. A linear theory shows that the transmitted wave angle is linked to the incident wave angle by a relation similar to Snell's law for surface waves (Pinettes et al., 1995). When the amplitude of the incident solitary wave is small, so that their celerity is very close to the critical phase speed, then we verify that the measured incident and transmitted wave angles are in agreement with the angles computed from Snell's law, using the critical phase speeds in the deep and shallow parts of the tank. However, when the solitary wave amplitude is large and the celerity significantly different from the critical phase speed, then the transmitted angle is smaller than predicted by the linear theory. The discrepancy between the predicted and the measured transmitted angle increases with increasing wave amplitude. We check that Snell's law is still valid granted that the (linear) critical phase speeds are replaced by the actual (nonlinear) celerities measured over the deep and shallow waters respectively. Then it appears that nonlinearity accounts for most of the discrepancy between the linear approximation and the measurements (fig. 4). This result is new, to our knowledge. Rotation has an influence upon the amplitude of the incident wave, which decreases when rotation increases, and upon the waves' celerity which tends to increase when rotation increases. We check that, with or without rotation, there is a maximum of wave amplitude, located over the slope. It correponds to the local perturbation which is necessary to link the incident and the transmitted waves.

BJ Solitary free surface wave interaction with high frequency free surface waves.

The experiments are conducted in the 36 m long glass flume equiped with 2 different wave makers. At one end solitary wave are generated by a piston wave maker with highly efficient law of motions prescribed by computer program.

The high frequency waves (frequencies between 1 Hz and 4 Hz, amplitudes approximately of 0. 5 cm) are produced by a wedge wave maker driven sinusoidally by an electrical motor. The interactions are either strong (both wave trains travelling in the same direction) or weak (both wave trains travelling in the opposite direction). The main phenomenon is the phase shifting between the short wave train before interaction and the wave train after interaction.

The phase shift produced by the interaction between a surface solitary wave and a short surface wave can be analyzed using very simple arguments. We recall that a fluid particle affected by a solitary wave travelling at constant speed is dispaced horizontally by a given length l whcih therefore induces net mass flux associated to such waves.

During the interaction with a short wave both motion combine. As a first approximation this combination may be regarded as linear. In consequence the short surface wave profile is shifted. The phase shift t.rp is the measure of this displacement between the wave train before and the wave train after the passage of the solitary wave. An approximate expression is thus:

t.rp = � kho IT v3 vhri
where h0 is the water depth, k the short surface wave number and ° A solitary wave amplitude. Note that Arp is independant of the short wave frequency. Nevertheless experiments suggest a slight increasing departure. from this law with decreasing frequency (see [START_REF] Morse | Contract MAS3-CT95-0027[END_REF].

CJ Interfacial solitary wave interaction with high frequency free surface waves.

This type of interaction is physically very close to the ones encountered in Ocean. Oceanic long nonlinear internal waves of characteristic length L of about 1 km are advantageously reproduced with solitary waves. On the other hand oceanic swells of wavelength A between 200 and 4 00 m are in our set-up produced as high frequency waves. The ratio of horizontal wave scales A/ L is of order 0. 5 in both the experiments and the ocean. For this topic the experiments are also conducted in the glass flume of LEGI. However 15 m of the total length was isolated and filled with brine and fresh water on the top. This produces an interface of finite thickness (of order 2cm

for a fractionnal jump in density of I%). The experiments are instrumented with two kinds of probes. Ultrasonic probes measure the interface displacement and capacitive probes measure free surface displacement. The wave generators are on one hand the wedge wave maker described above and an internal wave generator of the type used by Michallet & Barthelemy (1997, Exp. in Fluids, 22, pp 3 80-3 86 ). The wave interaction is measured with an array of 6 probes of each type. Both strong interactions (both wave trains travelling in the same directions) and weak interactions (waves travelling in the opposite directions) are performed. Experimental records consist of free surface and interface displacements against time at 6 locations in the flume 20 cm apart. For this set of experiments phase shifts were also expected. The weak velocity field associated to internal solitary waves compared with surface ones produces phase shifts too small to be measured with sufficient accuracy. However surface wave modulations are easily noticed. Strong interactions show that the surface wave is modulated with a decrease in amplitude as it interacts with the internal wave during strong interactions while the contrary is measured for weak interactions. This seems to agree with wave refraction on slowly varying currents theory as reported by Mei (1992, The Applied Dynamics of Ocean Surface Waves, World Scinetfic, p107). This experimental study should be supplemented in a near future by analytical modelling based on refraction-like 20 theories.

D) Numerical hydrodynamical modelling

The simulation of radar signatures of oceanic internal waves requires the knowledge of several physical mechanisms affecting the radar signal modulation. Among these are physical processes related to the internal wave dynamics itself, to its interaction with the ocean surface, and to the interaction between the sea surface and the electromagnetic radar waves, which constitutes the radar imaging of the ocean surface. At ItM an attempt is being made to describe all these physical processes by a set of different numerical models (the third one will be consideredin the last section): 1) A state-of-the-art internal wave model has recently been developed at IfM. This model, which has been applied to the simulation of the tidally induced dynamics in the straits of Gibraltar and Messina in the Mediterranean Sea, is capable of describing the generation and propagation of internal bores and their disintegration into trains of internal solitary waves [START_REF] Brandt | Internal waves in the Strait of Messina studied by a numerical model and synthetic aperture radar images from the ERS 1/2 satellites[END_REF]. Moreover, it allows for the numerical computation of realistic surface current fields associated with the internal wave dynamics.

2) The variable surface current fields calculated by the internal wave model are used as input data for the wave current interaction model of IfM. The wave current interaction model consists of a fast algorithm for the integration of the action balance equation according to weak hydrodynamic interaction theory. The existing software allows the employment of different source functions and parameterizations of the equilibrium waveheight spectrum and the relaxation rate. It can handle fully two-dimensional surface current fields as well as time-dependent surface currents. Also events like blocking and trapping of wave packets at certain positions or propagation of waves through shear currents can be handled without numerical problems.

E) SWARM

SWARM is the theoretical/modelling wing of the overall precommercial software and is based on past OCS experience and IfM's state-of-the-art knowledge of radar signature formation. (its scope goes beyond the hydrodynamical processes dealt with in this part and includes radar response). Its purpose is to generate numerical models representing ocean surface wave spectra at a number of surface spatial points throughout a straining current system. At the moment only 1-0 steady-state current systems are allowed (satisfactory for the surface distribution of current in internal tides, internal wave trains, swell waves, flow over bathymetry etc). However, the user has control over the distribution of wavevectors in the basic wave spectrum, and the distribution of spatial points. The physics of wave refraction is used to determine wave paths, with the dispersion relation for deep water gravity/capillary waves being employed. The energy content of each wavevector is determined by the physics of wave action, with a variety of wind input and wave breaking terms to be made available (for example Romeiser et al., 1997). Because of the number of parameters involved, it is possible to generate many wave spectra models for a given current distribution, e.g., just consider the multiplicity of models that may be generated using different wind speeds and wind directions. The number of radar signatures that may be predicted (radar wavelength, incidence angle, polarisation, imaging mechanism) becomes prodigious. A database of computed signatures, called a 'features data base', is therefore set up using the same cataloguing system as WASP.

IV. LABORATORY SAR DEVELOPMENT

A) Radar construction

The aim of this part of the MORSE project is to design, construct and test a high resolution W-band polarimetric imaging radar. The radar is designed for a resolution of approximately 6 cm in range, requiring a sweep bandwidth of 3 GHz, and 6 mm in azimuth. The maximum range is ver y small (7 m) and so the transmit power requirement is just I 0 mW with antenna gains of 15 dBi. The VCO sweep bandwidth is 3 GHz, the pulse repetition frequency 6 2 5 Hz and the image is formed in 1 second, during which time it is assumed that there is little movement within the imaged area. Antennas of 1 5 dBi gain should be suitable for these laboratory tank measurements, giving a SIN ratio of 20 dB or better across the image. Because of the very short range requirement, the IF frequency range is very low (1 OOs of kHz) and so may be sampled directly.

The FMCW technique is preferred to a pulsed radar for such high resolution at relatively short ranges. The design of the radar is simplified by using separate transmit and receive antennas, and this has the additional benefit of providing a polarimetric capability supporting HH, VY and HY modes though, of course, these are obtained with successive passes of the radar rather than simultaneously. To obtain azimuthal resolution of a similar order to the range resolution requires the use of aperture synthesis techniques with the radar moving along a rail. Usually the motion would be linear, but in this case a circular path is taken which somewhat complicates the azimuth processing. Aperture synthesis has been shown to be compatible with the FMCW technique (Griffiths,198 8) and a practical radar of this kind has in fact already been built and demonstrated by UCL at lower frequencies.

The beamwidth of the transmit and receive antennas should be sufficiently broad to provide full coverage of the swath in the vertical plane. In the horizontal plane the maximum achievable azimuth resolution is half the azimuth dimension of the antenna. For this application, horn-type antennas are used with apertures of approximately 6 mm x 4 mm and with beamwidths of approximately 35 ° and gains of approximately 1 5 dBi. This gives a maximum achievable azimuth resolution of 3 mm, though the particular geometry used means that this resolution is not realised, however an azimuth resolution superior to the range resolution is easy to achieve.

The radar velocity is 1 m/s requiring a minimum theoretical PRF of 333 Hz though again because of the unusual geometry a somewhat higher value is required to suppress the aliased doppler beams. VCO linearity is important, and is managed by standard calibration techniques in which the VCO frequency-against-voltage characteristic is measured and the required voltage-against-time characteristic derived and stored in digital form in RAM. Motion errors, which become significant at 0.8 mm (a quarter of a wavelength), are likely to be a more challenging problem, especially in, field trials where autofocus processing will certainly be required. It is planned to use corner reflectors in order to determine the exact path of the radar and correctly focus the aperture. A summary of the principal radar parameters is given in table l.

BJ First trials

During the first 1 2 months of the contract the radar design was established and the unit was successfully constructed at UCL. The VCO linearity was measured and an EPROM was programmed with the necessary calibration at ENSTBr. Initial measurements at UCL confirmed that the radar was working as expected although there appeared to be excessive mixer noise. Subsequent trials during the next 6 months at LEGI revealed that sensitivity was indeed below expectations due to the low-frequency mixer noise which is a fundamental weakness of the homodyne radar achitecture employed in the design. Subsequently, efforts have been made to reduce the mixer noise and its impact by increasing the PRF (by a factor of 4 ) which increases the IF frequencies and means that the mixer operates in a region which suffers less from low-frequency noise. The increased PRF results in reduced SIN ratio from the azimuth processing, but this can be restored simply by coherently averaging 4 adjacent looks (or, indeed more than 4 if a further improvement is required at the expense of azimuth resolution). These improvements are currently being finalised, after which it is planned to undertake further trials in both the stationary flume tank and the Coriolis turntable at LEGI.

V. RADAR BACKSCATTER MODELLING

The objective of this part is on one hand to further develop models for radar signatures generation, and on the other hand to address the inverse problem. This results in basic research as well as in software developpement.

A) Imaging theory

The output of the model described in section 3 . 4 , a two dimensional modulated wave spectrum for each grid point, is fed into the advanced composite surface scattering model of ItM for the numerical computation of radar signatures. This model is based on the standard expression of Bragg scattering theory, but it includes a number of second-order terms of the normalized radar backscattering cross section (NRCS) of the sea surface, which are obtained from a Taylor expansion of the NRCS with respect to the surface slopes parallel and normal to the radar look direction. The second-order terms of the NRCS account for the contribution of waves which are long compared to the short Bragg waves (Romeiser and Alpers, 1997;Romeiser et al., 1997). For each radar image according to the composite surface scattering model, an additional image is generated for comparison, which represents the result of a simple Bragg model without higher-order terms.

For the comparison between simulated and measured radar signatures of internal waves we used ERS 112 SAR images acquired over the straits of Gibraltar and Messina. It turned out that only very few of the available SAR images of internal waves are suited for a complete quantitative comparison with model results. This results from the fact that the signatures of internal waves are often noisy and two-dimensional in a way which is not accounted for in the internal wave model, or they are camouflaged by superimposed signatures of fronts, atmospheric features, or surface films. While the locations of the internal solitary waves does not seem to be affected strongly by such phenomena, the modulation depth of the NRCS can be affected very much. The analysis of internal wave signatures on ERS SAR images acquired over the Strait of Messina shows that the location of the internal solitary waves can be used to infer characteristics of the interior ocean. In fact, [START_REF] Brandt | Internal waves in the Strait of Messina studied by a numerical model and synthetic aperture radar images from the ERS 1/2 satellites[END_REF] show by using the internal wave model of IfM that the distance between the first two internal solitary waves of a wave train can be related, at least qualitatively, to the depth of the thermocline. Here, we want to focus on the possibility to gain information about the interior ocean from an analysis of the measured modulation depth of the NRCS over internal waves.

For first simulation runs an ERS 2 SAR image was selected, which shows signatures of an internal wave train near the Strait of Gibraltar. The signatures are well visible and have a regular shape. A comparison of the simulated surface current field with the observed radar signatures shows reasonable agreement between the measured and simulated locations of internal solitary waves. The wind speed at the station on the sh�re at the time of acquisition of the image was specified to be 6.1 m/s with a direction towards -1 5 ° with respect to our x-axis (which denotes the direction of the current as well as the azimuthal radar look direction). Fig. 6a shows the corresponding simulated NRCS profile (marked as 'l '), together with the measured one as obtained from the ERS 2 SAR image (marked as 'ERS'). In order to demonstrate that the simulated radar signatures depend significantly on the choice of the wind speed and direction, the figure includes also results for a wind direction rotated by 90° (line '2') and for a reduced wind speed of 3.0 mis (line '3'). In order to demonstrate the effect of the second order terms of the NRCS in the composite surface scattering model, Fig. 6b shows model results which would be obtained from a simple Bragg model without higher-order terms. ------------------, �8: �1 - Although the simulated and measured variations of the NRCS agree quite well without any dedicated tuning, the absolute (mean) level is clearly overestimated by the model in case 'I'. This behavior is already known from other simulated ERS SAR images and may result from the relatively steep incidence angle of 23 °, which is close to the limits of models based on Bragg scattering theory. However, our results for various wind speed vectors show that the absolute NRCS as well as the modulation depth predicted by the model depends clearly on the wind speed and direction, which is often not well known. We conclude from this that, even if we had a perfect imaging model, not only the complex structures in typical radar images but also the uncertainties of the available wind information would make it difficult to derive information on the surface currents of internal waves just from a comparison between measured and simulated modulation depths of the NRCS.

"I\ I � g-�- \ --- • 3 � 0 ���� 35 ----� 4 0 ��� ---'---' 45 Distance Fr om Sill (E.ast) [km] 0.8 �------------------, 0.7 0.6 � 0.

B) RADSIG

RADSIG combines WASP and SWARM and introduces a third component, MATCH. From WASP the user obtains an 'observed' signature (transect) from a desired and processed image; from S WARM the user has a data base of computed signatures made according to his/her intelligent input specification. By specifying/limiting parameters, the user enables MATCH to search the data base for the 'best fit' model that most closely resembles the observed signature. The 'inversion' process is the statement that the parameters of the observed signature are then those of the best-matched computed model. This inversion procedure has been followed based on our past experience and our knowledge of the complexity of the physics involved, e.g., internal hydrodynamics coupled with surface hydrodynamics and environmental conditions, coupled with electrodynamical backscattering from a near-random rough surface, and SAR technology. In brief, there is NO simple inversion technique and one is forced to the model-matching route. The algorithm presently intended for use in MATCH is based on matching the Fourier content of the observed signatures with the Fourier content of the computed signatures. This has the advantage that the Fourier data for the computed signatures may be computed but once and stored along side that of the 'physical' data. Preliminary studies of this physically motivated approach suggest that the method is quick, captures most of the essential information, and also restricts the number of 'matching' data base signatures to either one or a few, where the few may then be interrogated by the user. Most usefully, a user may 'customize' the feature data base so that it is representative of conditions that may be met with on a regional and seasonal basis, e.g., Straits of Gibraltar in September say.

C) Internal wave field reconstruction methodology

This part of the project is directly devoted to addressing the inverse problem of retrieving information on the water column from the manifestation of internal waves on radar images and is handled at TMS. This is a difficult inverse problem since highly non-linear transfer function are involved. More precisely, the corresponding forward problem consists in predicting radar signatures from some parameters describing the volume of the ocean. This forward problem decomposes into three successsive ones: The first one is the prediction of the surface current associated with the internal wave. The governing parameter is the density profile, especially the picnocline deph and density contrast. Assuming some external excitation, namely the tidal wave, the transfer function is given by a propagation model for the internal wave. The second one consists in evaluating the perturbed wave spectra for a given (straining) surface current, using weak hydrodynamic interaction theory. The first stage of the SW ARM software is dedicated to this task. The third one is the calculation of radar backscattering cross sections for a given surface wave spectrum. This can be done with various models ranging from simple Bragg scattering to advanced composite surface model as developped by IfM (see section 5.1 and 3.5).

One can search to solve the global inverse problem "directly" using a global forward model. It seems more realistic to adress successively the above mentionned problems, or at least two distinct ones, namely the first one and the other two together . This is made possible since they are independent in the sense that the output of one problem is the input of the following one. There are some advantages to proceed in this way. It is likely to save computer time with regards to the global method and allows to keep physical insight on the different physical processes involved. On the other hand, the method presents the drawback of relaxing some constraints on the solutions. In fact this is not a drawback if the model constraining the solution is not so reliable. This way to proceed has then been chosen for first investigation of the inverse problem. The main characteristic/difficulty of the inverse problem associated to the first forward problem is the multiplicity of models describing internal waves propagation either linear (modes, rays) or non-linear (Korteweg-De Vries, Benjamin Ono, Joseph), some of wich only qualitatively match the experimental data. This leads to a non standardly posed inverse problem. We have focussed on non-linear waves propagating far from their generation place. The main feature of the second forward problem is the large amount of computer time required to perform the calculation of perturbed surface wave spectra, especially when two-dimensionnal spectrum a r e taken into account. This precludes the use of global optimization schemes suitable for hihly non-linear inversion. At least as a first step, we limit our investigations to one dimensionnal spectra. The third forward problem, taken alone, would lead to an inverse problem that needs, to be solved, some minimization constraint. The development of works related to the inverse problem will be an important issue of the third year of the project.

VI. CONCLUSIONS

Substancial progress have been achieved towards the MORSE objective to better understand the physical processes involved in radar signatures of oceanic internal waves: On the experimental ground, nonlinear waves interaction experiments have been conducted and interpreted. For the refraction of internal solitary wave by a step shelf, a Snell's law is valid granted that nonlinear celerities are taken into account. In the case of long surface wave -short surface wave interaction, a phase shift of the short wave has been measured and shown to behave like predicted by an analytical model. Regarding SAR image processing, besides classical algorithms for parameters extraction, an automatic method for internal wave detection and location using wavelet transform has been successfully applied to images from two important areas for internal waves generation. The project dedicated high resolution W-band radar has been designed, constructed and tested. First measurements show that the radar is working as expected but with an excessive mixer noise. Improvements are being finalised. On the numerical modelling ground, an improved model for computing two-dimensional perturbed wave spectra has been developped. Once associated to a composite surface model for backcatter, it leads to the computation of radar signatures. The influence on these signatures of wind speed and direction, often badly known parameters, has been emphasised. Three user-friendly precommercial software packages are produced during the project, version 1 of the first two being completed: WASP for extraction of wave parameters from SAR images, SWARM for radar signatures generation, and RADSIG for matching observed signatures with computed ones.

1

  TMS: Thomson Marconi Sonar, 7-9 rue des Mathurin$, 9 2221 Bagneux cedex, France 2 IfM: Institut fur Meereskunde, Universitat Hamburg 3 LEG!: Laboratoire des Ecoulements geophysiques et industriels 4 UCL: University College London 5 OCS: Oxford Computer Services ltd 6 ENSTB: Ecole Nationale Superieure des Telecommunications de Brest • Development of methods for the extraction of ocean surface characteristics at high spatial resolution and for the resolution of the inverse problem consisting of retrieving the internal wave fields from surface information.

  are presented in the corresponding sections. They constitute the elements of a pre-competitive global tool for automatic extraction of information on ocean hydrological and dynamical conditions from radar images. As part of the data management plan of the project a MORSE Web site ll has been set up at OCS. This web site has been building up gradually and now contains general information regarding MORSE activities, some ERS-1 imagery of internal waves, and samples of data from the experiments carried out in the flume tank and Coriolis tank by LEGI scientists. Interested users may download and analyse some of this sample data from the web site before requesting complete data sets from LEGI (please visit the web site for further information).

Fig. 1 :

 1 Fig. 1: Stages of the internal wave detection and location system several packets of internal solitons together with radar signatures of possible tropical rain cells and dark spots due to wind conditions. The fact of having several islands close to the coastline introduces a problem in the detection system. Fig. 3(b) show the results obtained from the application of the internal wave detection method. Most of the internal wave features have been detected although the difficulty due to noise, dark spots, and tropical rain cells in the image. It can be observed that the detected internal wave crests have their origin in the interaction with the shallow bottom topography of the shelf, which coincides with the estimations in the Andaman Sea given in (Alpers et al., 1997). 3) Perspectives An automatic method for internal wave detection and location has been applied to two satellite SAR images, and the feasibility of this technique for feature extraction has been demonstrated. The presented coastline detection method is accurate enough for the application on the internal wave detection problem. The wavelet transfonn in its multiscale edge detection version has been seen to be an excellent method for internal wave feature extraction. Several SAR images from two of the main world's places in generation of internal waves in the sea surface have been studied. Internal wave structures are well detected and located against background noise by means of local modulus maxima methods across scales. Internal wave packets are quite well distinguished against internal waveJook-alikes using an a priori oceanographic knowledge of the internal wave geometry shape and size.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Internal wave feature detection. (a) ERS-1 SAR image from the Strait of Gibraltar which was acquired on September IO, 1993, covering an area of 100 km x JOO km. It shows several internal wave packets together with other oceanic features. (b) Detected internal waves features by the proposed method.

Fig. 4 :

 4 Fig. 4: Refraction angle as function of the incident angle. Open symbols (o and 0) correspond to the measured angles for two different initial conditions and the x• and + to the angles computed from the actual nonlinear celerities.

  5i:s::: : 0. 4 .

Fig. 6 :

 6 Fig. 6: Simulated surface current (bottom) and radar signatures for a scenario corresponding to an existing ERS 2 SAR image of internal waves near the Strait of Gibraltar. 'ERS' -NRCS as obtained from the ERS 2 SAR image, '1' -NRCS as obtained from the model for the nominal wind speed (6.1 m/s) and direction (-15° with respect to the x axis), '2' -NRCS as obtained from the model for a wind direction rotated by 90°, '3' -NRCS as obtained from the model for a reduced wind speed of 3.0 mis; (a) full composite surface scattering model, (b) Bragg model without higher-order terms. The azimuthal radar look direction is equal to the x-direction.

Table 1 .

 1 Principal radar parameters

	radar centre frequency: 9 4 GHz sweep bandwidth 3 GHz PRF 6 2 5 Hz imaging time 1 s transmit power l OmW radar velocity 1 m/s size of antennas 6 mm x 4 mm gain of antennas 1 5 dBi range resolution 6 cm azimuth resolution 1. 5 cm
	operating range	7m
	SIN ratio	>20 dB

ACKNOWLEDGEMENTS

The Morse consortium gratefully acknowledges the CEC MAST Program for its support for this project under contract MAS3-CT-95-0027 (DG XII). The ENSTB group would like to thank the Centre for Remote Sensing and Processing, CRISP, National University of Singapore, for providing a WWW Catalogue Browse whence to obtain the RASARSA T SAR images from the Andaman Sea.

The research is being undertaken by a Consortium of four University and two commercial partners under the MAST