Vincent Giovangigli
email: giovangi@cmap.polytechnique.fr

Milan Pokorný
email: pokorny@karlin.mff.cuni.cz

Ewelina Zatorska
email: e.zatorska@mimuw.edu.pl

On the steady flow of reactive gaseous mixture

Keywords: multicomponent flow, chemically reacting gas, steady compressible Navier-Stokes-Fourier system, weak solution

We present the study of systems of equations governing a steady flow of polyatomic, heat-conducting reactive gas mixture. It is shown that the corresponding system of PDEs admits a weak solution and renormalized solution to the continuity equation, provided the adiabatic exponent for the mixture γ is greater than 5 3 .

Introduction

In this paper we consider a system of equations modelling the steady flow of compressible, heat-conducting, n-component mixtures, undergoing reversible chemical reactions in a bounded domain Ω ⊂ R 3 . These are n reaction-diffusion equations coupled to the Navier-Stokes-Fourier system describing the fluid motion. The applications of such systems are numerous, especially in engineering, like spatial vehicle reentry, crystal and polymers growth, combustion or atmospheric pollution [START_REF] Giovangigli | Multicomponent flow modeling[END_REF], but they may also be used to describe population or the chemotaxis models [START_REF] Chen | Analysis of a multidimensional parabolic population model with strong cross-diffusion[END_REF][START_REF] Gajewski | Global behaviour of a reaction-diffusion system modelling chemotaxis[END_REF]. This is a strong motivation for investigating the mathematical structure and properties of the corresponding PDEs. Our goal is to extend the existence theory [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF][START_REF] Zatorska | On a steady flow of multicomponent, compressible, chemically reacting gas[END_REF] to the case of heat-conducting mixtures with strong cross-diffusion.

The motion of gaseous mixture can be described by a system of equations governing the total mass density = (x), the velocity vector field u = u(x), the absolute temperature ϑ = ϑ(x) and the species concentrations Y k (x), k ∈ {1, . . . , n}. These equations express the conservation of the total mass, momentum, total energy and species masses. They may be written as div( u) = 0, div( u ⊗ u) + divS + ∇π = f , div( Eu) + div(πu) + divQ + div(Su) = f • u, div( Y k u) + divF k = m k ω k , k ∈ {1, . . . , n}.

(

In the above equations S denotes the viscous tensor (it has the opposite sign to the viscous stress tensor also considered in continuum mechanics), π the internal pressure of the fluid, f the external force, E the specific total energy, Q the heat flux, ω k the molar production rate of the k-th species, F k the diffusion flux of the k-th species and m k the molar mass of the k-th species.

The fundamental difficulty in such type of systems is the coupling between the fluid mechanics part, governed by the Navier-Stokes-Fourier system, with the species reactiondiffusion equations. The coupling appears through the form of the pressure, which may depend on the species concentrations, and through the heat-flux which, in contrast to the single-component flows, includes a term representing the transfer of energy due to species molecular diffusion.

Here we focus on the second type of coupling, while we assume that the molar masses are comparable, i.e.

m 1 = . . . = m n = 1, (2)
which leads to the pressure independent of the species concentrations. Assumption [START_REF] Bruno | Relaxation of quantum state population and volume viscosity in he/h2 mixtures[END_REF] is fulfilled for the mixtures of isomers or in the state-to-state fluid models where each quantum state is a separate pseudo species [START_REF] Bruno | Relaxation of quantum state population and volume viscosity in he/h2 mixtures[END_REF]. Similar assumption was made in [START_REF] Feireisl | Multicomponent reactive flows: global-in-time existence for large data[END_REF],

where the existence of variational weak entropy solutions was shown, yet for the diagonal Fick diffusion. This approximation, however, does not take into account the cross-effects that are well-known to play an important role in the flows of multicomponent fluids. In the recent papers [START_REF] Zatorska | On a steady flow of multicomponent, compressible, chemically reacting gas[END_REF][START_REF] Zatorska | Multicomponent mixture model. the issue of existence via time discretization[END_REF] the model of isothermal flow with the Fick diffusion and pressure depending on the species concentration was considered. However, it turns out that this assumption leads to inconsistency with the Second Law of Thermodynamics when the heat-conductivity is taken into account.

Reaction-diffusion equations with more complex diffusion and even their coupling to the incompressible Navier-Stokes equations were investigated using many different approaches [START_REF] Bothe | On the Maxwell-Stefan approach to multicomponent diffusion[END_REF][START_REF] Bulíček | On existence of weak solution to a model describing incompressible mixtures with thermal diffusion cross effects[END_REF][START_REF] Chen | Analysis of an incompressible Navier-Stokes-Maxwell-Stefan system[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF][START_REF] Herberg | Reaction-diffusion systems of maxwell-stefan type with reversible mass-action kinetics[END_REF][START_REF] Marion | Global existence for fully nonlinear reaction-diffusion systems describing multicomponent reactive flows[END_REF]. Our goal in this paper is to provide analogue existence result with no restriction on the size of data but in the compressible setting. So far, such results were available only for system with viscosity coefficients satisfying the so called BD relation and some further restrictions on the form of the pressure [START_REF] Mucha | Heat-conducting, compressible mixtures with multicomponent diffusion: construction of a weak solution[END_REF][START_REF] Zatorska | On the flow of chemically reacting gaseous mixture[END_REF].

In our case, due to relatively weak coupling between the Navier-Stokes equations and the rest of system (1), a part of techniques is a combination of methods developed first by Lions [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF] and Feireisl [START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF] for evolutionary barotropic flows and by Mucha, Pokorný [START_REF] Mucha | On the steady compressible Navier-Stokes-Fourier system[END_REF] and Novotný, Pokorný [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] for the stationary Navier-Stokes-Fourier system. The main difference and the biggest difficulty in the present paper concerns much more complex form of entropy inequality, which is a source of majority of a-priori estimates. Therefore, the approximation scheme had to be considerably modified in comparison to the single-fluid case considered in [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF]. Some ideas for construction of the solution to the subsystem of species reaction-diffusion equations can be already found in our previous work [START_REF] Mucha | Chemically reacting mixtures in terms of degenerated parabolic setting[END_REF]. Here, however, we had to extend them by adding further regularizations necessary to handle general form of diffusion matrix.

The outline of the paper is the following. In the next section we introduce the main assumptions on the parameters of our model. In Section 3 we present the definition of the weak solution and formulate the main result: existence of a weak solution under certain restriction on the parameters of our model. Sections 4 and 5 contain the description of the approximate scheme. In Section 6 we prove existence of a solution to the full approximation and in the subsequent sections, 7 and 8, we pass to the limit with all regularizing parameters to get existence of a solution to the original problem.

Formulation of the problem

We consider system (1) supplemented by the no-slip boundary conditions

u| ∂Ω = 0, (3)
together with

F k • n| ∂Ω = 0, (4)
and the Robin boundary condition for the heat flux

-Q • n + L(ϑ -ϑ 0 ) = 0. (5)
The last condition means that the heat flux through the boundary is proportional to the difference of the temperature inside Ω and the known external temperature ϑ 0 . We assume that the total mass of the mixture is given

Ω dx = M > 0. ( 6
)
The mass fractions Y k , k ∈ {1, . . . , n} are defined by Y k = k , where k are the species densities and n k=1 k = . Thus, by definition, the mass fractions satisfy

n k=1 Y k = 1. (7)
The consistency with the principle of mass conservation requires that the diffusion fluxes and the species production rates satisfy

n k=1 F k = 0, n k=1 ω k = 0, (8)
keeping in mind that the species molar masses are equal (2).

Fundamental thermodynamic relations

We consider a pressure π = π( , ϑ) of the following form

π( , ϑ) = π c ( ) + π m ( , ϑ), (9)
where π m obeys the Boyle law

π m = n k=1 Y k ϑ = ϑ, (10)
i.e., it represents the pressure for an ideal mixture of n species, whose molar masses are equal 1 and we take, without loss of generality, the gaseous constant R = 1. The first component of ( 9), π c , is the so called cold pressure or the barotropic correction

π c = γ , γ > 1.
The specific total energy E is a sum of the specific kinetic and specific internal energies

E( , u, ϑ, 1 , . . . , n ) = 1 2 |u| 2 + e( , ϑ, Y 1 , . . . , Y n ),
where the latter, similarly to the pressure, consists of two components

e = e c ( ) + e m (ϑ, Y 1 , . . . , Y n ).
The cold energy e c and the ideal gas mixture energy e m are given by

e c = 1 γ -1 γ-1 , e m = n k=1 Y k e k = ϑ n k=1 c vk Y k ,
where c vk is the mass constant-volume specific heat. The constant-pressure specific heat, denoted by c pk , is related (under assumption (2)) to c vk in the following way

c pk = c vk + 1, (11)
and both c vk and c pk are assumed to be constant (but possibly different for each constituent).

In agreement with the second law of thermodynamics, there exists a differentiable function called the specific entropy of the mixture s( , ϑ, Y 1 , . . . , Y n ) that can be expressed in terms of the partial specific entropies

s k = s k ( , ϑ, Y k ) of the k-th species s = n k=1 Y k s k . ( 12
)
The differential of entropy is related to the differential of energy, total density and mass fractions by the Gibbs formula

ϑDs = De + πD 1 - n k=1 g k DY k , (13)
with the Gibbs functions

g k = h k -ϑs k , (14)
where

h k = h k (ϑ), s k = s k ( , ϑ, Y k
) denote the specific enthalpy and the specific entropy of the k-th species, respectively, with the following exact forms

h k = c pk ϑ, s k = c vk log ϑ -log -log Y k .
The cold pressure and the cold energy correspond to an isentropic processes, therefore using (13) one can derive an equation for the specific entropy s

div( su) + div Q ϑ - n k=1 g k ϑ F k = σ, ( 15
)
where σ is the entropy production rate

σ = -S : ∇u ϑ - Q • ∇ϑ ϑ 2 - n k=1 F k • ∇ g k ϑ - n k=1 g k ω k ϑ . ( 16
)

The form of transport fluxes

The viscous tensor S is determined by the Newton rheological law as

S = S(ϑ, ∇u) = -µ ∇u + (∇u) t - 2 3 divuI -ν(divu)I, (17)
where µ = µ(ϑ) > 0, ν = ν(ϑ) ≥ 0 are the shear and bulk viscosity coefficients, respectively, I is the identity matrix.

The heat flux Q consists of two terms (we already neglect the Soret and the Dufour effects) representing the transfer of energy due to the species molecular diffusion and the Fourier law, respectively

Q = n k=1 h k F k + q, q = -κ∇ϑ, (18)
where κ = κ(ϑ) > 0 is the thermal conductivity coefficient. The diffusion flux of the k-th species F k is given by

F k = -Y k n l=1 D kl ∇Y l , (19)
where

D kl = D kl ( , ϑ, Y 1 , . . . , Y n ), k, l = 1, . . . , n are the multicomponent diffusion coeffi- cients.
The diffusion matrix. The coefficients D kl depend only on ϑ and Y 1 , . . . , Y n , see for instance [START_REF] Giovangigli | Multicomponent flow modeling[END_REF], therefore we introduce another matrix

(D kl ) n k,l=1 = (D kl ) n k,l=1 = (D kl (ϑ, Y 1 , . . . , Y n )) n k,l=1 .
The main properties of the diffusion matrix D, discussed in [START_REF] Giovangigli | Multicomponent flow modeling[END_REF], Chapter 7, are:

D = D t , N (D) = R Y , R(D) = Y ⊥ , D is positive semidefinite over R n , (20)
where we assumed that Y = (Y 1 , . . . , Y n ) t > 0. Above N (D) denotes the nullspace of matrix D, R(D) denotes its range, U = (1, . . . , 1) t and U ⊥ denotes the orthogonal complement of RU . Furthermore, we assume that the matrix D is homogeneous of a non-negative order with respect to Y 1 , . . . , Y n and that D ij are differentiable functions of ϑ, Y 1 , . . . , Y n for any i, j ∈ {1, . . . , n} such that

|D ij (ϑ, Y )| ≤ C( Y )(1 + ϑ a )
for some a ≥ 0.

Remark 1 As a consequence of (20) the matrix D is positive definite over U ⊥ . This property corresponds to the positivity of entropy production rate associated with the diffusive process, see [START_REF] Waldmann | Formale kinetische Theorie von Gasgemischen aus anregbaren Molekülen[END_REF]. Indeed, according to above definitions, σ may be rewritten in the following form σ = -S(ϑ, ∇u) :

∇u ϑ + κ|∇ϑ| 2 ϑ 2 - n k=1 F k • ∇ (log p k ) - n k=1 g k ω k ϑ , (21)
where we denoted p k = Y k ϑ. Let us investigate the structure of the third term, we have

- n k=1 F k • ∇ (log p k ) = - n k=1 F k p k • ∇p k = - n k=1 F k • ∇Y k Y k + ∇( ϑ) ϑ [due to (8)] = n l,k=1 D kl ∇Y l • ∇Y k ≥ c n k=1 |∇Y k | 2 . ( 22
)
The transport coefficients. The coefficients µ, ν, κ are continuous functions of temperature and the following growth conditions are imposed

µ(1 + ϑ) ≤ µ(ϑ) ≤ µ(1 + ϑ), 0 ≤ ν(ϑ) ≤ ν(1 + ϑ), κ(1 + ϑ m ) ≤ κ(ϑ) ≤ κ(1 + ϑ m ), (23)
with m > 0 and the positive constants µ, µ, ν, κ, κ.

The species production rates. We assume that the species production rates are smooth bounded functions of ( , ϑ, Y 1 , . . . , Y n ) such that

ω k ( , ϑ, Y 1 , . . . , Y n ) ≥ 0 whenever Y k = 0. (24)
Next, in accordance with the second law of thermodynamics we assume that

- n k=1 g k ω k ≥ 0, ( 25
)
where g k are specified in [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF]. Note that thanks to this inequality, [START_REF] Novotný | Weak solutions for steady compressible Navier-Stokes-Fourier system in two space dimensions[END_REF], together with [START_REF] Mucha | Chemically reacting mixtures in terms of degenerated parabolic setting[END_REF] and [START_REF] Mucha | Approximate solutions to model of two-component reactive flow[END_REF] yield that the entropy production rate defined in ( 16) is non-negative.

Notation

The matrices from R n×n are denoted by capital letters, the vectors from R n are denoted by . The vectors from R 3 and the tensors from R 3×3 are denoted by small and capital bold letters, respectively. We use generic constant denoted by C which may change from line to line. When it is important, its dependence of parameters will be indicated in the parentheses.

We work in the framework of Sobolev and Lebesgue spaces denoted by W m,p (Ω), m ∈ N, and L p (Ω), p ≥ 1, respectively, endowed with the standard norms. For brevity we will write

u W m,p (Ω) = u m,p , u L p (Ω) = u p , u L p (∂Ω) = u p,∂Ω
independently whether u is a vector or scalar. By C ∞ 0 (Ω) we denote the space of C ∞ functions on Ω with zero value at the boundary ∂Ω.

Weak solutions, main result

We are now in a position to formulate the definition of weak solutions for our system. Definition 1 We say the set of functions ( , u, ϑ, Y ) is a weak solution to problem (1-5), [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF][START_REF] Feireisl | Multicomponent reactive flows: global-in-time existence for large data[END_REF][START_REF] Gajewski | Global behaviour of a reaction-diffusion system modelling chemotaxis[END_REF][START_REF] Giovangigli | Multicomponent flow modeling[END_REF][START_REF] Herberg | Reaction-diffusion systems of maxwell-stefan type with reversible mass-action kinetics[END_REF][START_REF] Lions | Mathematical topics in fluid mechanics[END_REF][START_REF] Marion | Global existence for fully nonlinear reaction-diffusion systems describing multicomponent reactive flows[END_REF][START_REF] Mucha | On the steady compressible Navier-Stokes-Fourier system[END_REF][START_REF] Mucha | Chemically reacting mixtures in terms of degenerated parabolic setting[END_REF][START_REF] Mucha | Approximate solutions to model of two-component reactive flow[END_REF], [START_REF] Waldmann | Formale kinetische Theorie von Gasgemischen aus anregbaren Molekülen[END_REF][START_REF] Zatorska | On a steady flow of multicomponent, compressible, chemically reacting gas[END_REF]

provided ≥ 0 a.e. in Ω, ∈ L γ (Ω), Ω dx = M , u ∈ W 1,2 0 (Ω), |u| 3 ∈ L 1 (Ω), Y k ≥ 0 a.e. in Ω, Y ∈ W 1,2 (Ω), F k • n| ∂Ω = 0, and n k=1 Y k = 1 a.e. in Ω, ϑ > 0, a.e. in Ω, ϑ m ∇ϑ ∈ L 1 (Ω), ϑ ∈ L 1 (∂Ω),
- Ω (u ⊗ u) : ∇ϕ ϕ ϕ + S : ∇ϕ ϕ ϕ dx - Ω πdivϕ ϕ ϕ dx = Ω f • ϕ ϕ ϕ dx, holds for any test function ϕ ϕ ϕ ∈ C ∞ 0 (Ω); • the weak formulation of the species equations - Ω Y k u • ∇ψ dx - Ω F k • ∇ψ dx = Ω ω k ψ dx,
holds for any test function ψ ∈ C ∞ (Ω) and for all k = 1, . . . , n;

• the weak formulation of the total energy balance

- Ω 1 2 |u| 2 + e u • ∇ψ dx + Ω κ∇ϑ • ∇ψ dx - Ω n k=1 h k F k • ∇ψ dx = Ω f • uψ dx + Ω (Su) • ∇ψ dx + Ω πu • ∇ψ dx - ∂Ω L(ϑ -ϑ 0 )ψ dS, holds for any test function ψ ∈ C ∞ (Ω).
Remark 2 Indeed, there exists a more general definition of a solution to our problem which is based on replacing the weak formulation of the total energy balance by the entropy inequality and the global total energy balance; see [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] or [START_REF] Novotný | Weak and variational solutions to steady equations for compressible heat conducting fluids[END_REF] for more details. The reason for this is that it allows to prove the existence of those variational entropy solutions for larger range of parameters. Nonetheless, we prefer to stay here with (a more direct) definition of weak solutions to our problem, even for the prize of more restrictive conditions on m and γ.

Remark 3

The kinetic energy balance is (at least formally) nothing but the momentum equation multiplied by u and integrated over Ω, thus we may write the balance of the internal energy in the form

div( eu) + div n k=1 h k F k -κ∇ϑ + πdivu + S : ∇u = 0. ( 26
)
However, the balance of the total energy is equivalent to the balance of the internal energy and the momentum equation only for sufficiently regular solutions, but it might be not true for weak solutions introduced above.

We will also use the notion of the renormalized solution to the continuity equation

Definition 2 Let u ∈ W 1,2 loc (Ω) and ∈ L 6/5 loc (Ω) solve div( u) = 0
in the sense of distributions on Ω, then the pair ( , u) is called a renormalized solution to the continuity equation, if

div(b( )u) + ( b ( ) -b( ))divu dx = 0, ( 27
)
in the sense of distributions on Ω, for all b ∈ W 1,∞ (0, ∞) ∩ C 1 ([0, ∞)), such that sb (s) ∈ L ∞ (0, ∞).
The main theorem of this paper reads as follows

Theorem 1 Let γ > 5 3 , M > 0, m > 1, a < 3m-2 2 . Let Ω ∈ C 2 .
Then there exists at least one weak solution to our problem above in the sense of Definition 1. Moreover, ( , u) is the renormalized solution to the continuity equation in the sense of Definition 2.

First level of approximation

This section is devoted to the main level of approximation on which all essential approximation parameters appear: ε > 0 indicating additional dissipation and relaxation in the continuity equation as well as in the species mass balance equations, λ > 0 providing elliptic regularization in the species equations while written in terms of entropic variables, and δ > 0 improving the integrability of the density and providing the bound from below for the temperature. Combining ideas from [START_REF] Feireisl | Multicomponent reactive flows: global-in-time existence for large data[END_REF][START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF][START_REF] Zatorska | On a steady flow of multicomponent, compressible, chemically reacting gas[END_REF] we consider the following approximative system. We look for ∈ W 2,q (Ω) ∩ L β (Ω) for some q

≥ 6 5 , u ∈ W 1,2 0 (Ω), Y ∈ W 1,2 (Ω), ϑ ∈ W 1,2 (Ω) ∩ L 3B (Ω) ∩ L 2B (∂Ω) with ≥ 0, Y k ≥ 0 and ϑ > 0 a.e.
in Ω satisfying the following system.

• The approximate continuity equation

ε + div( u) = ε∆ + ε , ∇ • n| ∂Ω = 0, (28)
is satisfied pointwisely and we require > 0, Ω dx = M , thus we may take = M |Ω| . • The weak formulation of the approximate momentum equation

Ω 1 2 u • ∇u • φ φ φ - 1 2 (u ⊗ u) : ∇φ φ φ -S : ∇φ φ φ dx - Ω (π + δ β + δ 2 )divφ φ φ dx = Ω f • φ φ φ dx (29) is satisfied for each φ φ φ ∈ C ∞ 0 (Ω).
• The weak formulation of the approximate species balance equations (k = 1, 2, . . . , n)

Ω εY k ψ dx - Ω Y k u • ∇ψ dx - Ω F k • ∇ψ dx + λ Ω (∇ log Y k • ∇ψ + log Y k ψ) dx = Ω ω k ψ -ε ∇Y k • ∇ψ + εdiv(Y k ∇ )ψ -ε∇Y k • ∇ψ + ε k ψ dx (30)
is satisfied for any ψ ∈ C ∞ (Ω), where k > 0, k = 1, . . . , n satisfy n k=1 k = ; for instance, we take k = n = M n|Ω| . Moreover,

F k = -Y k n l=1 D kl (ϑ, Y )∇Y l , D kl (ϑ, Y ) = 1 (σ Y + ε) r D kl (ϑ, Y ), (31)
for suitable r ≥ 1 (connected with the order of the homogeneity of

D(•, Y )) and σ Y = n k= Y k .
• The weak formulation of the approximate total energy equation

- Ω e + 1 2 |u| 2 + (π + δ β + δ 2 ) u • ∇ψ dx - Ω Su • ∇ψ + δϑ -1 ψ dx + Ω κ δ ε + ϑ ϑ ∇ϑ • ∇ψ dx + ∂Ω (L + δϑ B-1 )(ϑ -ϑ 0 ) + ε log ϑ + λϑ B 2 log ϑ ψ dS + n k=1 c vk Ω -F k • ∇ψ + ϑ(ε( + 1)Y k + λ) ∇Y k Y k • ∇ψ dx = Ω f • uψ dx + δ β -1 Ω (εβ β-1 ψ + β u • ∇ψ -εβ β ψ) dx + δ Ω (2ε ψ + 2 u • ∇ψ -εβ 2 ψ) dx (32)
is satisfied for any ψ ∈ C ∞ (Ω), where

κ δ = κ + δϑ B + δϑ -1 .
Remark 4 Above β and B are some positive, large enough numbers that will be determined in the course of the proof.

For this system we will prove the following result Theorem 2 Let 0 < λ ε δ, β and B be sufficiently large positive numbers. Let Ω ∈ C 2 . Then there exists a solution to the approximate system in the sense specified above.

The existence of solutions to the above system will be proven below, by introducing another artificial level of approximation.

Full approximation

Now, our task is to construct regular solution defined in the previous section, i.e. to prove Theorem 2. For this purpose we introduce two new parameters: N ∈ N -the dimension of the Galerkin approximation in the momentum equation and η -regularization of the coefficients in the temperature and the momentum equations. Note that since the species molar masses are assumed to be equal, the species concentrations do not appear in the momentum equation. This allows to solve more or less separately the Navier-Stokes-Fourier system and the reaction-diffusion system and to combine them via a suitable fixed point theorem. In order to consider the former, we used slightly modified strategy from [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] by Novotný and Pokorný. The existence of weak solutions to a similar system of reactiondiffusion equations in the evolutionary case is due to Mucha, Pokorný and Zatorska [START_REF] Mucha | Chemically reacting mixtures in terms of degenerated parabolic setting[END_REF]. In [START_REF] Mucha | Approximate solutions to model of two-component reactive flow[END_REF], [START_REF] Mucha | Heat-conducting, compressible mixtures with multicomponent diffusion: construction of a weak solution[END_REF] the same authors investigated also the coupling between the two systems, note however, that the assumptions on the form of the fluxes and pressure are now quite different. Nevertheless we use some of their arguments here.

In the basic level of approximation we look for the set of functions { , u, Y , ϑ} satisfying the following system.

• The approximate continuity equation

ε + div( u) = ε∆ + ε , ∇ • n| ∂Ω = 0 (33)
is satisfied pointwisely.

• The Galerkin approximation for the momentum equation

Ω 1 2 u • ∇u • w - 1 2 (u ⊗ u) : ∇w -S η : ∇w dx - Ω (π + δ β + δ 2 )divw dx = Ω f • w dx (34) is satisfied for each test function w ∈ X N , X N = span{w i } N i=1 ⊂ W 1,2 0 (Ω), where {w i } N i=1
are the first N eigenfunctions of the Laplace operator with Dirichlet boundary condition; indeed, u ∈ X N . In (34) we denoted

S η = - µ η (ϑ) 1 + ηϑ ∇u + (∇u) t - 2 3 divuI - ν η (ϑ) 1 + ηϑ (divu) I. (35)
• The approximate species mass balance equations

divJ k = R k , J k • n| ∂Ω = 0 (36)
are satisfied pointwisely. In (36) we denoted

J k = - n l=1 Y k Y l D kl (ϑ, Y )∇Y l /Y l -ε( + 1)Y k + λ ∇Y k /Y k , R k = ω k ( , ϑ, Y ) + ε k -εY k -div(Y k u) + εdiv(Y k ∇ ) -λ log Y k , (37)
and D kl (ϑ, Y ) was defined in (31).

• The approximate internal energy balance

-div κ δ,η ε + ϑ ϑ ∇ϑ = -div( eu) -πdivu + δ ϑ -S η : ∇u + δε(β β-2 + 2)|∇ | 2 -div ϑ n k=1 c vk J k , (38)
with the boundary condition

κ δ,η ε + ϑ ϑ ∇ϑ • n| ∂Ω + (L + δϑ B-1 )(ϑ -ϑ η 0 ) + ε log ϑ + λϑ B 2 log ϑ = 0, (39)
is satisfied pointwisely and ϑ η 0 stands for a smooth, strictly positive approximation of ϑ 0 .

Remark 5 Note that in fact we expect to have in (38) in the last term c pk instead of c vk , which can be easily observed by setting λ = ε = δ = 0. However, since n k=1 F k = n k=1 F k = 0, both c vk and c pk lead after the limit passages to the same. Hence we prefer to keep c vk in the approximate scheme.

In the above system κ δ,η , µ η , ν η are regularizations of functions κ δ , µ, and ν extended by constants κ δ (0), µ(0) and ν(0) to the negative half-line.

The existence of solutions is formulated in the following theorem Theorem 3 Let δ, ε, λ and η be positive numbers and N be a positive integer. Let Ω ∈ C 2 . Then there exists a solution to system (33-36) such that

∈ W 2,q (Ω), ∀q < ∞, ≥ 0 in Ω, Ω dx = M , u ∈ X N , Y ∈ W 1,2 (Ω) with log Y k ∈ W 2,q (Ω) ∀q < ∞, Y k ≥ 0 a.e. in Ω and ϑ ∈ W 2,q (Ω), ∀q < ∞, ϑ ≥ C(N ) > 0.
The strategy of the proof is the following:

1. we rewrite the system above for τ = log ϑ and

Z k = log Y k , k = 1, 2, . . . , n;
2. we fix u in the space X N and use it to find a unique smooth solution to (33) = (u); here we may follow [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF] verbatim;

3. we find a unique solution to a system of equations which is nothing but a linearization of the above system (34)-(39) written in the new variables;

4. we apply a Schauder type of fixed point theorem for the momentum, the internal energy and the species equations and we deduce the existence of u ∈ X N and log ϑ ∈ W 2,q (Ω) log Y k ∈ W 2,q (Ω); this part follows similarly to [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] provided some a-priori estimates are valid.

Existence of solutions for the full approximation

Step 1: We define the operator

S : X N → W 2,p (Ω), 1 ≤ p < ∞, S(u) =
, where solves the approximate continuity equation (28) with the Neumann boundary condition. We then claim Lemma 4 Let assumptions of Theorem 3 be satisfied. Then the operator S is well defined for all p < ∞. Moreover, if S(u) = , then ≥ 0 in Ω and

Ω dx = Ω dx = M . Additionally, if u X N ≤ L, L > 0, then 2,p ≤ C(ε, p, Ω, M )(1 + L), 1 < p < ∞. ( 40
)
The above lemma is an analogue of Proposition 4.29 from [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF], so we omit the proof.

Step 2: We consider the mapping:

T : X N × W 1,∞ (Ω) × W 1,∞ (Ω) → X N × W 1,∞ (Ω) × W 1,∞ (Ω) with ( S = (S 1 , . . . , S n ) t , Z = (Z 1 , . . . , Z n ) t ) T (v, S, σ) = (u, Z, τ ),
where (we take always = (v) the unique solution to (33) with u replaced by v)

Ω -S η (e τ , ∇u) : ∇w dx = Ω 1 2 (v ⊗ v) : ∇w - 1 2 v • ∇v • w dx + Ω π( , e σ ) + δ β + δ 2 divw dx + Ω f • w dx, (41)
which is fulfilled for any w ∈ X N ,

-div n l=1 D kl (e σ , e S )e S k e S l ∇Z l + ε( + 1)e S k + λ ∇Z k + λZ k =ω k + ε -ε e S k -div( ve S k ) + εdiv(e S k ∇ ), k = 1, . . . , n (42)
) 44
The existence of a unique solution to (41-44) is a consequence of the Lax-Milgram theorem. Since the r.h. sides and the boundary terms are of lower order and sufficiently smooth, the operator T is compact. The continuity of T is straightforward.

In what follows we will show that there exists a bound (independent of t) for all fixed points to

tT (u, Z, τ ) = (u, Z, τ ) (45)
for t ∈ [0, 1]. Thus, due to a well known version of the Schauder fixed point theorem (see f.i. [7, Theorem 9.2.4]) we finish the proof of Theorem 3. In fact, we shall show that these estimates are independent of N , η and λ which will be used in the subsequent limit passages.

Uniform estimates

We now denote ϑ = e τ , thus ϑ > 0, and we can in particular divide the energy equation by ϑ. Similarly we denote Y k = e Z k and have also Y k > 0, but we do not know yet if n k=1 Y k = 1. However, we denote

σ Y = n k=1 Y k
and in what follows, we show certain estimates of σ Y -1 (we get the norm small provided λ is small). Take t ∈ [0, 1] and consider the fixed points as in (45), i.e. consider:

• the continuity equation

ε + div( u) = ε∆ + ε , ∇ • n| ∂Ω = 0, (46)
• the momentum equation

Ω -S η : ∇w dx = t Ω 1 2 (u ⊗ u) : ∇w - 1 2 u • ∇u • w dx + t Ω (π + δ β + δ 2 )divw dx + t Ω f • w dx, (47)
• the species mass balance equations

-div n l=1 Y k D kl ∇Y l + ε( + 1)Y k + λ ∇Y k Y k + λ log Y k = t ω k + ε k -ε Y k -div( uY k ) + εdiv(Y k ∇ ) , k = 1, . . . , n, (48)
• the internal energy equation

-div κ δ,η ε + ϑ ϑ ∇ϑ = -tdiv( eu) -tπdivu + tδε(β β-2 + 2)|∇ | 2 -tS η : ∇u + t δ ϑ + n k=1 c vk div ϑ n l=1 Y k D kl ∇Y l + ϑ ε( + 1)Y k + λ ∇Y k Y k , (49)
with the boundary conditions

n l=1 Y k D kl ∇Y l + ε( + 1)Y k + λ ∇Y k Y k • n| ∂Ω = 0 (50) κ η,δ ε + ϑ ϑ ∇ϑ • n| ∂Ω + t(L + δϑ B-1 )(ϑ -ϑ η 0 ) + ε log ϑ + λϑ B 2 log ϑ = 0. ( 51
)
We first add up the equations for the species (48) for k = 1, . . . , n and test the obtained sum by ξ ∈ C ∞ (Ω), we get

ε Ω ( + 1)∇σ Y • ∇ξ dx + λ n k=1 Ω ∇ log Y k • ∇ξ dx + λ n k=1 Ω log Y k ξ dx = t Ω ε -ε σ Y -div( uσ Y ) + εdiv(σ Y ∇ ) ξ dx.
On the other hand, the continuity equation can be written as

ε Ω ξ dx - Ω u • ∇ξ dx = -ε Ω ∇ • ∇ξ dx + ε Ω ξ dx.
Therefore, multiplying the latter by t and subtracting the two equations above we may write

t Ω ε(σ Y -1) ξ -(σ Y -1)u • ∇ξ + ε(σ Y -1)∇ ∇ξ dx + ε Ω ( + 1)∇σ Y • ∇ξ dx = -λ n k=1 Ω ∇ log Y k • ∇ξ dx -λ n k=1 Ω log Y k ξ dx. Now, taking ξ = (σ Y -1) we obtain ε Ω ( + 1)|∇σ Y | 2 dx - Ω u • ∇(σ Y -1) 2 2 dx + t Ω ε (σ Y -1) 2 + ε(σ Y -1)∇ • ∇(σ Y -1) dx = -λ n k=1 Ω ∇ log Y k • ∇σ Y dx -λ n k=1 Ω log Y k (σ Y -1) dx. ( 52
)
Multiplying the continuity equation by

(σ Y -1) 2 2 we obtain ε (σ Y -1) 2 2 + div( u) (σ Y -1) 2 2 = ε∆ (σ Y -1) 2 2 + ε (σ Y -1) 2 2 , (53)
thus, multiplying (53) by t, integrating it over Ω and subtracting from (52), we see that

tε Ω (σ Y -1) 2 2 dx + tε Ω (σ Y -1) 2 2 dx + ε Ω ( + 1)|∇(σ Y -1)| 2 dx = -λ n k=1 Ω ∇ log Y k • ∇σ Y dx -λ n k=1 Ω log Y k (σ Y -1) dx. ( 54
)
Hence we get an estimate for n

k=1 Y k ε t σ Y -1 2 2 + ∇σ Y 2 2 ≤ C(ε)λ 2 n k=1 ∇ log Y k 2 2 + 1 t log Y k 2 2 , (55)
which will become very important in derivation of estimates following from the entropy inequality.

Next, using w = u in (47) we obtain

- Ω S η : ∇u dx = t Ω [π + δ( β + 2 )]divu + f • u dx. ( 56
)
Integrating (49) over Ω and using (51), we get that

∂Ω t(L + δϑ B-1 )(ϑ -ϑ η 0 ) + ε log ϑ + λϑ B 2 log ϑ dS = -t Ω S η : ∇u dx + t Ω δϑ -1 dx -t Ω πdivu dx + tεδ Ω |∇ | 2 (β β-2 + 2) dx.
We add up these two expressions and use the renormalized continuity equation, to get

∂Ω t(L + δϑ B-1 )(ϑ -ϑ 0 ) + ε log ϑ + λϑ B 2 log ϑ dS + tεδ Ω β β -1 β + 2 2 dx -(1 -t) Ω S η : ∇u dx + ε ∂Ω log ϑ dS =t Ω f • u dx + t Ω δ ϑ dx + tεδ Ω β β -1 β-1 + 2 dx. (57)
Moreover, dividing (49) by ϑ we derive

-div κ δ,η ε + ϑ ϑ 2 ∇ϑ -κ δ,η (ε + ϑ) ϑ |∇ϑ| 2 ϑ 2 = -t u ∇e ϑ - π ϑ ∇ 2 -tdiv( u) e ϑ + π ϑ + n k=1 c vk div n l=1 Y k D kl ∇Y l + (ε + 1)Y k + λ ∇Y k Y k - n k,l=1 c vk Y D kl ∇Y l • ∇ log ϑ - n k=1 c vk ε( + 1)Y k + λ ∇Y k Y k • ∇ log ϑ -t S η : ∇u ϑ + t δ ϑ 2 + tδε(β β-2 + 2) |∇ | 2 ϑ . ( 58
)
Now we multiply (48) by ξ = log Y k , sum up with respect to k and integrate over Ω

n k,l=1 Ω D kl ∇Y k • ∇Y l dx + n k=1 Ω ε( + 1) |∇Y k | 2 Y k + λ|∇ log Y k | 2 + λ(log Y k ) 2 dx =t n k=1 Ω ω k + ε k -ε Y k -div( uY k ) + εdiv(Y k ∇ ) log Y k dx. ( 59
)
Next we multiply (48) on ξ = c vk log ϑ -c pk , sum up with respect to k and integrate over Ω:

n k,l=1 c vk Ω Y k D kl ∇Y l • ∇ log ϑ dx + n k=1 c vk Ω ε( + 1)Y k + λ ∇Y k Y k • ∇ log ϑ dx = n k=1 λ Ω log Y k c pk -c vk log ϑ dx + t n k=1 c vk Ω ω k + ε k -ε Y k -div( uY k ) + εdiv(Y k ∇ ) log ϑ dx -t n k=1 c pk Ω ω k + ε k -ε Y k -div( uY k ) + εdiv(Y k ∇ ) dx. ( 60
)
Now, take (57)-Ω (58)dx-( 59)-(60):

n k,l=1 Ω D kl ∇Y k • ∇Y l dx + n k=1 Ω ε( + 1) |∇Y k | 2 Y k + λ|∇ log Y k | 2 dx + λ n k=1 Ω (log Y k ) 2 dx + Ω κ δ,η (ε + ϑ) ϑ |∇ϑ| 2 ϑ 2 dx + t Ω -S η : ∇u ϑ + δ ϑ 2 + δε(β β-2 + 2) |∇ | 2 ϑ dx -(1 -t) Ω S η : ∇u dx + tεδ Ω β β -1 β + 2 2 dx + t ∂Ω (L + δϑ B-1 )ϑ dS + t ∂Ω (L + δϑ B-1 ) ϑ η 0 ϑ dS + ∂Ω∩{ϑ>1} ε log ϑ + λϑ B 2 log ϑ dS - ∂Ω∩{ϑ≤1} ε log ϑ ϑ + λϑ B 2 -1 log ϑ dS = RHS,
where

RHS = t ∂Ω (L + δϑ B-1 )ϑ η 0 dS + t ∂Ω (L + δϑ B-1 ) dS - ∂Ω∩{ϑ≤1} ε + λϑ B 2 ϑ dS + ∂Ω∩{ϑ>1} εϑ -1 + λϑ B 2 -1 log ϑ dS + t Ω f • u dx + t Ω δ ϑ dx + tεδ Ω β β -1 β-1 dx + tεδ Ω 2 dx + t Ω u ∇e ϑ - π ϑ ∇ 2 dx + t Ω div( u) e ϑ + π ϑ dx + t n k=1 Ω ω k log Y k + ε k log Y k -ε Y k log Y k -div( uY k ) log Y k dx + t n k=1 Ω εdiv(Y k ∇ ) log Y k dx -t n k=1 c vk Ω ω k log ϑ + ε k log ϑ dx + t n k=1 c vk Ω ε Y k log ϑ + div( uY k ) log ϑ -εdiv(Y k ∇ ) log ϑ dx + λ n k=1 c vk Ω log Y k log ϑ dx -λ n k=1 c pk Ω log Y k dx + t n k=1 c pk Ω ω k + ε k -ε Y k dx = 25 i=1 I i . (61)
Finally, from (56) and the renormalized continuity equation we see

u 2 1,2 + tεδ( β β + ∇ β/2 2 2 ) + tεδ( 2 2 + ∇ 2 2 ) ≤ t Ω (πdivu + f • u) dx
and from the standard elliptic regularity for the continuity equation we also have

∇ q ≤ C(ε)(1 + |u| q ), 1 < q < ∞.
Hence,

u 2 1,2 + tεδ( β β + ∇ β/2 2 2 ) + tεδ( 2 2 + ∇ 2 2 ) ≤ Ct 2 (1 + ϑ 2β β-2 2β β-2 ), ( 62
)
provided β > 2, and

∇ q ≤ C(ε) 1 + t ϑ β β-2 2β β-2 , 1 < q ≤ 6β β + 2 . ( 63
)
We can now estimate the terms on the r.h.s. of (61). First, note that we may easily estimate the boundary integrals and the terms I 5 -I 8 . In order to estimate terms I 9 and I 10 , recall that we may write

e = 1 γ -1 γ-1 + ϑ n k=1 c vk Y k = e c + e m .
We first estimate the part corresponding to e c . We have

t Ω u ∇e c ϑ - π ϑ ∇ 2 dx + t Ω div( u) e c ϑ + π ϑ dx = -t Ω u • ∇ dx + t Ω div( u) 1 γ -1 γ-1 ϑ + 1 + γ-1 ϑ dx = t Ω div( u) log + γ γ -1 γ-1 ϑ dx = t Ω (ε∆ + ε -ε ) log + γ γ -1 γ-1 ϑ dx = -εt Ω |∇ | 2 dx -εt { ≥1} log dx + εt { <1} log dx -εt { <1} log dx + εt { ≥1} log dx -εt Ω γ ϑ γ γ -1 dx -tεγ Ω γ-2 |∇ | 2 ϑ dx + εt γ γ -1 Ω γ-1 ∇ • ∇ϑ dx + εt γ γ -1 Ω γ-1 ϑ dx.
Now, observe that the terms no. 1, 2, 3, 6, 7 from the r.h.s. of above are non-positive, therefore they can be transferred to the l.h.s. of (61). The terms no. 4, 5, 8, 9 can be then easily estimated by means of those terms and the rest of the terms from the l.h.s. of (61). This part of the proof is a repetition of arguments from [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] so we skip the details.

Next we take I 9 (e m ) + I 10 (e m ) + I 14 + I 19 and get

t n k=1 Ω u ϑ • ∇(ϑc vk Y k ) + div( u)(c vk Y k ) -div( uY k ) log Y k dx + n k=1 c vk Ω log ϑdiv( uY k ) dx = t Ω u • ∇σ Y dx = -t Ω (ε∆ -ε + ε )σ Y dx.
Then, to handle the first term we use ( 55) and ( 62)

εt Ω ∇ • ∇σ Y dx ≤ εt ∇ 2 ∇σ Y 2 ≤ C(ε)t 1 + ϑ β β-2 2β β-2 n k=1 λ ∇ log Y k 2 + 1/ √ t log Y k 2 ,
which can be estimated by the l.h.s. of (61). The second term can be estimated similarly, and the last term is non-positive. Now, for I 11 + I 16 + I 23 , recalling that n k=1 ω k = 0, we obtain

t n k=1 Ω ω k (log Y k -c vk log ϑ + c pk ) dx =t n k=1 Ω ω k (log Y k -c vk log ϑ + c pk + log ) dx = t Ω n k=1 g k ω k ϑ dx ≤ 0.
Further, for I 12 we use again (55) to write

I 12 ≤ t n k=1 Ω∩{Y k ≤1} ε k log Y k dx + Ω∩{Y k >1} ε k log Y k dx ≤ t n k=1 Ω∩{Y k >1} ε k Y k dx ≤ C 1 + λ n k=1 ( √ t ∇ log Y k 2 + log Y k 2 )
and the terms on the r.h.s. of above can be now easily estimated by the l.h.s. of (61). Next, I 13 may be split into two parts, namely

I 13 = -tε n k=1 Ω∩{Y k >1} Y k log Y k dx + Ω∩{Y k ≤1} Y k log Y k dx ≤ Ctε Ω dx ≤ C(M ).
The term I 15 can be estimated similarly as I 9 + I 10 for e m . The term I 17 can be easily estimated by means of t Ω δ ϑ 2 dx from the l.h.s. of (61). For α > 0 small, we may write

I 18 = εt Ω log ϑ n k=1 c vk Y k dx ≤ Cεt Ω | log ϑ|(σ Y -1) dx + Ctε Ω | log ϑ| dx ≤ Cεt β ( σ Y -1 2 + 1) ∂Ω ϑ B 2 log ϑ dS + Ω κ δ,η (ϑ) |∇ϑ| 2 ϑ 2 dx + ε ∂Ω∩{ϑ<1} log ϑ ϑ dS α ≤ λ 4 n k=1 ( ∇ log Y k 2 2 + log Y k 2 2 ) + λ 4 ∂Ω ϑ B 2 log ϑdS + ε 4 Ω κ δ,η (ϑ)|∇ϑ| 2 dx + ε 4 ∂Ω∩{ϑ<1} log ϑ ϑ dS.
To treat I 20 , we will use the continuity equation to replace the highest order terms of the density

I 20 = -tε n k=1 c vk Ω div(Y k ∇ ) log ϑ dx = -tε n k=1 c vk Ω ∇Y k • ∇ log ϑ dx -tε n k=1 c vk Ω Y k ∆ log ϑ dx ≤ Ctε Ω n k=1 |∇Y k | 2 Y k dx 1 2 σ Y -1 1 2 2 ∇ 4 log ϑ 6 -t n k=1 c vk Ω Y k (ε -ε + div( u)) log ϑ dx ≤ λ 4 n k=1 ( ∇ log Y k 2 2 + log Y k 2 2 ) + ε n k=1 ∇ Y k 2 2 + λ 4 ∂Ω ϑ B 2 log ϑ dS + ε 4 ∂Ω∩{ϑ<1} log ϑ ϑ dS + tε 4 ∂Ω ϑ B dS + ε 4 Ω κ δ,η |∇ϑ| 2 dx,
and

I 21 = λ n k=1 Ω c vk log Y k log ϑ dx ≤ λ 4 n k=1 log Y k 2 2 + λ 4 ∂Ω ϑ B 2 log ϑdS + ε 4 Ω κ δ,η (ϑ)|∇ϑ| 2 dx + ε 4 ∂Ω∩{ϑ<1}
log ϑ ϑ dS, so, I 18 , I 20 , I 21 may be estimated by the l.h.s. of (61). The terms I 22 , I 24 and I 25 are easy or can be estimated as above. Now we may employ the standard elliptic theory (cf. [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF]) to show that we can estimate log ϑ and log Y k in W 2,q (Ω) for any q < ∞, i.e. in particular, in W 1,∞ (Ω) independently of t. The proof of Theorem 3 is finished.

2 7 Proof of Theorem 2 7.1 Limit passage N → ∞
Recall that in the previous section, we proved the following estimates (see (61) with t = 1). Note that the constant on the r.h.s. is independent of N , η and λ, however, may depend on ε and δ:

√ λ Y 1,2 + n k=1 |∇Y k | 2 Y k 1 + ∇ϑ B 2 2 + ∇ϑ ϑ 2 2 + ϑ -2 1 + ϑ B,∂Ω + log ϑ ϑ 1,∂Ω + β ≤ C. (64)
Moreover, (62) and (63) together with (33) imply

∇ 2 2 + u 1,2 + ∇ 6 ≤ C. (65)
Hence, we easily pass with N → ∞ in the continuity equation (33), in the weak formulation of the momentum equation (34) and the weak formulation of the species mass balance balance equations (36). However, we cannot pass to the limit so easily in the internal balance equation due to the presence of the term S η (ϑ, ∇u) : ∇u which is bounded only in L 1 (Ω). Nonetheless, we may proceed as in [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF]: we may use as test function in the limit momentum balance the function u to show that lim

N →∞ Ω S η (ϑ N , ∇u N ) : ∇u N dx = Ω S η (ϑ, ∇u) : ∇u dx
which yields the strong convergence of the velocity gradient in L 2 (Ω) and the passage to the limit in the weak formulation of the internal energy balance (38) can be performed. Let us stress that this is the only limit passage for which the strong convergence of the velocity gradient can be shown. It is caused by the fact that the momentum equation may be tested by the solution u only when η > 0.

7.2 Limit passage η → 0 + Hence, before the next limit passage η → 0 + we must replace the internal energy balance by the total energy balance. Note that at this level, they are still equivalent -we may test the momentum equation by u -but this might not be true after the limit passage η → 0 + . To this aim, we use in the momentum equation as a test function uψ with ψ ∈ C ∞ (Ω) and sum it with the weak formulation of the internal energy balance with the test function ψ. We get the weak formulation of the approximate total energy balance

- Ω e + 1 2 |u| 2 + (π + δ β + δ 2 ) u • ∇ψ dx - Ω S η u • ∇ψ + δϑ -1 ψ dx + Ω κ δ,η ε + ϑ ϑ ∇ϑ • ∇ψ dx + ∂Ω (L + δϑ B-1 )(ϑ -ϑ η 0 ) + ε log ϑ + λϑ B 2 log ϑ ψ dS + n k=1 c vk Ω ϑ n l=1 Y k D kl ∇Y l • ∇ψ + ϑ(ε( + 1)Y k + λ) ∇Y k Y k ψ dx = Ω f • uψ dx + δ Ω ε|∇ | 2 (β β-2 + 2) + ( β + 2 )divu ψ dx.
Therefore, using the renormalized continuity equation we deduce

- Ω e + 1 2 |u| 2 + (π + δ β + δ 2 ) u • ∇ψ dx - Ω S η u • ∇ψ + δϑ -1 ψ dx + Ω κ δ,η ε + ϑ ϑ ∇ϑ • ∇ψ dx + ∂Ω (L + δϑ B-1 )(ϑ -ϑ η 0 ) + ε log ϑ + λϑ B 2 log ϑ ψ dS + n k c vk Ω ϑ n l=1 Y k D kl ∇Y l • ∇ψ + ϑ(ε( + 1)Y k + λ) ∇Y k Y k • ∇ψ dx = Ω f • uψ dx + δ β -1 Ω (εβ β-1 ψ + β u • ∇ψ -εβ β ψ) dx + δ Ω (2ε ψ + 2 u • ∇ψ -2ε 2 ψ) dx
for all ψ ∈ C ∞ (Ω). Hence we may now use ( 64)-(65) to justify the limit passage η → 0 + and hence to verify the validity of ( 28)-(32). Theorem 2 is proved. 2

Proof of Theorem 1

To finish the proof of the main theorem we subsequently pass with λ → 0 + , ε → 0 + and δ → 0 + .

8.1 Limit passage λ → 0 +

The first limit passage in this section is still based on estimates (64)-(65). Note, however, that we cannot dispose with the estimate of Y in W 1,2 (Ω) anymore. On the other hand, we still have from (55)

σ Y -1 6 ≤ C(λ), (66)
with C(λ) → 0 when λ → 0 + . Therefore, combining estimate (66) with the second term on the l.h.s. of (64) we get

∇ Y 12 7 ≤ C, (67)
with C independent of λ. Therefore, using (64), ( 65), ( 66) and (67) we may pass with λ → 0 + and get • the approximate continuity equation

ε + div( u) = ε∆ + ε , ∇ • n| ∂Ω = 0 (68)
• the weak formulation of the approximate momentum equation

Ω 1 2 u • ∇u • φ φ φ - 1 2 (u ⊗ u) : ∇φ φ φ -S : ∇φ φ φ dx - Ω (π + δ β + δ 2 )divφ φ φ dx = Ω f • φ φ φ dx, (69)
satisfied for each φ φ φ ∈ C ∞ 0 (Ω) • the weak formulation of the approximate species balance equations

ε Ω Y k ψ dx - Ω Y k u • ∇ψ dx + Ω n l=1 Y k D kl ∇Y l • ∇ψ dx = Ω ω k ψ -ε ∇Y k • ∇ψ + εdiv(Y k ∇ )ψ -ε∇Y k • ∇ψ + ε k ψ dx, (70) satisfied for any ψ ∈ C ∞ (Ω) (k = 1, 2, . . . , n)
• the weak formulation of the approximate total energy equation

- Ω e + 1 2 |u| 2 + (π + δ β + δ 2 ) u • ∇ψ dx - Ω Su • ∇ψ + δϑ -1 ψ dx + Ω κ δ ε + ϑ ϑ ∇ϑ • ∇ψ dx + ∂Ω (L + δϑ B-1 )(ϑ -ϑ 0 ) + ε log ϑ ψ dS + Ω ϑ n k,l=1 c vk Y k D kl ∇Y l • ∇ψ + ϑ n k=1 ε( + 1)c vk ∇Y k • ∇ψ dx = Ω f • uψ dx + δ β -1 Ω (εβ β-1 ψ + β u • ∇ψ -εβ β ψ) dx + δ Ω (2ε ψ + 2 u • ∇ψ -2ε 2 ψ) dx (71)
satisfied for any ψ ∈ C ∞ (Ω).

Limit passage ε → 0 +

At this step we loose all the control of the density, except for the L 1 -bound due to the given mass. On the other hand, due to (66), we know now that σ Y = 1 which leads to the bounds (recall

Y k ≥ 0) 0 ≤ Y k ≤ 1, ∀k ∈ {1, 2, . . . , n}.
Unlike the situation before, we may also deduce from the first term on the lhs of (61) additional bound on ∇ Y . Employing again σ Y = 1, we have

Ω n k,l=1 D kl ∇Y k ∇Y l dx ≥ C ∇ Y 2 2 .
Further, due to the Korn inequality and the form of the viscous tensor, we control the W 1,2 -norm of the velocity. We can therefore estimate most of the terms on the r.h.s. of (61) and get (cf. [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF])

Y 1,2 + Y ∞ + ∇ϑ B 2 2 + ϑ B,∂Ω + ϑ 3m + ϑ -2 1 + ϑ -1 1,∂Ω + u 1,2 ≤ C 1 + Ω f • u dx . (72)
However, due to lack of sufficient density estimate, the boundedness of the r.h.s. of (72) has to be verified. Note that the momentum equation is in fact the same as in the case of the compressible Navier-Stokes-Fourier system studied in [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF], so, we may apply the same technique to obtain the so called Bogovskii-type of estimates. Following [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF], we use as test function in (69) the function φ φ φ, solution to

divφ φ φ = 2 3 β - 1 |Ω| Ω 2 3 β dx, φ φ φ| ∂Ω = 0.
For more information on the Bogovskii operator, we refer the reader to [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF], Lemma 3.17 and to [START_REF] Danchin | The divergence equation in rough spaces[END_REF]. In consequence of this testing we may obtain the additional bound on , namely

- Ω ϑ + γ + δ( β + 2 ) divφ φ φ dx = Ω f • φ φ φ dx, (73)
for all φ φ φ ∈ C ∞ 0 (Ω) • the weak formulation of the approximate species balance equations

- Ω Y k u • ∇ψ dx + Ω Y k n l=1 D kl (ϑ, Y )∇Y l • ∇ψ dx = Ω ω k ψ dx, for all ψ ∈ C ∞ (Ω), k = 1 = 2, . . . , n
• the weak formulation of the approximate total energy equation

- Ω ϑ n k=1 c vk Y k + 1 2 |u| 2 + ( ϑ + γ γ -1 γ + δ( β + 2 )) u • ∇ψ dx - Ω Su • ∇ψ + δϑ -1 ψ dx + Ω κ δ ∇ϑ • ∇ψ dx + ∂Ω (L + δϑ B-1 )(ϑ -ϑ 0 )ψ dS + Ω ϑ n k,l=1 c vk Y k D kl ∇Y l • ∇ψ dx = Ω f • uψ dx + δ Ω 1 β -1 β + 2 u • ∇ψ) dx (74)
satisfied for any ψ ∈ C ∞ (Ω). Now, and in the sequel, the bar denotes the weak limit of the corresponding term, for instance γ = lim ε→0 + γ ε in the sense of the weak convergence in L 5β 3γ (Ω). Recall that since ∈ L 2 (Ω), we also have (see e.g. [START_REF] Feireisl | Singular limits in thermodynamics of viscous fluids[END_REF]) the continuity equation satisfied in the renormalized sense (cf. Definition (2)).

Next, we should also show the strong convergence of the density. Since this can be shown exactly as in [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF] and we meet similar problems in the next subsection when passing with δ → 0 + , we skip here all details. Hence we may remove the bars in (73) and (74).

Limit passage δ → 0 +

Last subsection is devoted to the limit passage δ → 0 + , i.e. to the final step of the proof of Theorem 1. Similarly as in the previous section, we may use the entropy inequality to deduce a-priori estimates independent of δ. However, in this case, we must proceed more carefully.

The main idea is that unlike the previous limit passage, we subtract the total energy balance tested by a constant function from our entropy estimate and deduce estimate of the type

∇ Y 2 2 + ∇ϑ m 2 2 2 + u 2 1,2 + ϑ -1 1,∂Ω + δ ∇ϑ B 2 2 2 + ∇ϑ -1 2 2 2 + ϑ -2 1 + ϑ B-2 1,∂Ω ≤ C(1 + δ ϑ B-1 1,∂Ω ). ( 75
)
Recall that as σ Y = 1, we also have

Y ∞ ≤ C.
From the total energy balance tested by a constant function we deduce

ϑ 1,∂Ω + δ ϑ B 1,∂Ω ≤ C 1 + Ω u • f dx + δ ϑ -1 1 . (76)
To get rid of the δ-dependent terms we use once more the Bogovskii-type of estimates, i.e. we test the momentum equation by φ φ φ -a solution to

divφ φ φ = - M |Ω| , φ φ φ| ∂Ω = 0.
This allows us to get a bound δ β+1 β+1 ≤ C which can be employed to get rid of the δ-dependent terms in (75) and (76). It yields

∇ Y 2 2 + Y ∞ + ∇ϑ m 2 2 2 + u 2 1,2 + ϑ -1 1,∂Ω + δ( ∇ϑ B 2 2 2 + ∇ϑ -1 2 ) 2 2 + ϑ -2 1 + ϑ B-2 1,∂Ω ≤ C (77)
and

ϑ 3m ≤ C(1 + 6 5
).

The details can be again found in [START_REF] Novotný | Steady compressible Navier-Stokes-Fourier system for monoatomic gas and its generalizations[END_REF]. Thus we need now additional δ-independent estimates of the density. To this aim, we employ again the Bogovskii-type estimates i.e. we test the momentum equation by φ φ φ -a solution to . Using these bounds we may pass to the limit in our system to get the weak formulation of the continuity equation Note that here we used that δ converges weakly to in L γ+α (Ω), ϑ δ converges strongly to ϑ in L q (Ω) for any q < 3m. Again, γ stands for the weak limit of γ δ in L γ+α γ (Ω). Note also that all the δ-dependent terms tend strongly in L 1 (Ω) to zero. It is also not so difficult to pass to the limit in the species balance equations to get

divφ φ φ = α - 1 |Ω| Ω α dx, φ φ φ| ∂Ω = 0. We obtain Ω π( , ϑ) + δ( β + 2 ) α dx = - Ω (u ⊗ u) : ∇φ φ φ dx - Ω S(ϑ, ∇u) : ∇φ φ φ -f • φ φ φ dx + 1 |Ω| Ω p( , ϑ) + δ( β + 2 ) dx
- Ω Y k u • ∇ψ dx + Ω Y k n l=1 D kl ∇Y l • ∇ψ dx = Ω ω k ψ dx,
for all ψ ∈ C ∞ (Ω), k = 1, 2, . . . , n, under the assumptions that the growth in the D kl term with respect to temperature is below 3m 2 . Anyway, we get a stronger restriction below. We get the weak formulation of the total energy balance in the form . In order to pass to the limit in the term δ |u δ | 2 u δ we must require that the density is bounded in L p (Ω) for some p > 2 and to pass to the limit in the term with the stress tensor we must require that ϑ δ → ϑ strongly in L q (Ω) for q > 3. This leads to the restrictions m > 1 and γ + α > 2, i.e. γ > 5 3 and m > 1. In order to finish the proof, we have to remove the bars over certain nonlinear quantities which requires that δ → strongly in L 1 (Ω).

The last step of the proof of main theorem is hence to show the strong convergence of the density. Here we follow the ideas due to P.-L. Lions [START_REF] Lions | Mathematical topics in fluid mechanics[END_REF] developed for the isentropic flows, used in the context of the heat conducting fluid in [START_REF] Mucha | On the steady compressible Navier-Stokes-Fourier system[END_REF] for constant viscosities and in [START_REF] Novotný | Weak solutions for steady compressible Navier-Stokes-Fourier system in two space dimensions[END_REF] for temperature dependent viscosities. Due to the uniform L 2 bounds of the density, we can directly employ the renormalized solution to the continuity equation and no truncations in the test functions are needed.

We first verify the validity of the effective viscous flux identity in the form ( ϑ + γ ) α -(µ(ϑ) + ν(ϑ)) α divu = ϑ + γ α -(µ(ϑ) + ν(ϑ)) α divu, with α > 0 from the estimates above. It can be shown exactly as in Lemma 8 in [START_REF] Novotný | Weak and variational solutions to steady equations for compressible heat conducting fluids[END_REF] (with α instead of Θ). It is based on testing the momentum equation before and after limit passage by ξ(x)∇∆ -1 (1 Ω α δ ) and by ξ(x)∇∆ -1 (1 Ω α ), respectively, and on proving certain limit passages via compensated compactness technique. Next we may verify (note that both and δ belong to L 2 (Ω)) that the following versions of the renormalized continuity equation hold true div( α δ u δ ) + (α -1) α δ divu δ = 0 and div( α u δ ) + (α -1) α divu = 0, both in the sense of distributions in R 3 . Hence (see [START_REF] Novotný | Introduction to the mathematical theory of compressible flow[END_REF]Lemma 4.39]) div(( α )

1 α u) = 1 -α α ( α ) 1 α -1 ( α divu -α divu).
Therefore div((( α )

1 α -)u) = 1 -α α ( α ) 1 α -1 ( ϑ + γ ) α -ϑ + γ α µ(ϑ) + ν(ϑ) , (79)
again in the sense of distributions in R 3 . Testing (79) by ψ = 1 reads

Ω ( α ) 1 α -1 ( ϑ + γ ) α -ϑ + γ α µ(ϑ) + ν(ϑ) dx = 0.
It is easy to verify that δ → 0 in L 1 ({ α = 0}). Thus the monotonicity of t → t γ and the strong convergence of the temperature yield γ+α + ϑ 1+α = γ α + ϑ α .

As ϑ > 0 a.e. in Ω, we obtain α = α which implies 1+α = 1+α a.e. in Ω.

Whence the strong convergence of the density. The proof of Theorem 1 is now complete. 2

  and the following integral equalities hold • the weak formulation of the continuity equation Ω u • ∇ψ dx = 0, holds for any test function ψ ∈ C ∞ (Ω); • the weak formulation of the momentum equation

and- 2 + n k=1 c vk div e σ n l=1 e

 2l=1 div (κ δ,η (e σ ) (ε + e σ ) ∇τ ) = -div e( , e σ , e S )v -π( , e σ )divv-S η (e σ , ∇v) : ∇v + δe -σ + δε(β β-2 + 2)|∇ | S k e S l D kl (e σ , e S )∇Z l + e σ ε( + 1)e S k + λ ∇Z k ,with the boundary conditions n l=1 D kl (e σ , e S )e S k e S l ∇Z l + (ε( + 1)e S k + λ)∇Z k • n| ∂Ω = 0 (43) κ η,δ (e σ ) (ε + e σ ) ∇τ • n| ∂Ω + (L + δe σ(B-1) )(e σ -ϑ η 0 ) + ετ + λe qσ τ = 0. (

5 3 β

 3 ≤ C. Thus we may pass to the limit in (68)-(71) and get • the weak formulation of the continuity equation Ω u • ∇ψ dx = 0, for all ψ ∈ C ∞ (Ω) • the weak formulation of the approximate momentum equation -Ω ( (u ⊗ u) : ∇φ φ φ + S : ∇φ φ φ) dx

3 , 2 ≤ 2 ,

 322 the estimates of the most restrictive terms. The convective term yields Ω (u ⊗ u) : ∇φ φ φ dx ≤ C γ+α u which leads to the restriction 0 < α ≤ 2γ -3 for γ > 3 2 . The stress tensor can be estimated Ω S(ϑ, ∇u) : ∇φ φ φ dx ≤ C(1 + ϑ 3m ) u 1,2 φ φ φ 6m 3mwhich leads to restriction 0 < α ≤ 3m-2 3m+2 γ, m > 2 3 . Thanks to this we get, in addition to (77) and (78) γ+α ≤ C, where 0 < α ≤ min 2γ -3, 3m -2 3m + 2 γ with γ > 3 2 and m > 2 3

Ωu

  • ∇ψ dx = 0, for all ψ ∈ C ∞ (Ω), and in the weak formulation of the approximate momentum equation-Ω ( (u ⊗ u) : ∇φ φ φ + S : ∇φ φ φ) dx -Ω ϑ + γ divφ φ φ dx = Ω f • φ φ φ dx,for all φ φ φ ∈ C ∞ 0 (Ω).

c

  vk Y k D kl ∇Y l • ∇ψ dx = Ω f • uψ dx satisfied for any ψ ∈ C ∞ (Ω), under several additional restrictions. First of all, we must assume that D kl (ϑ, •) ≤ C(1 + ϑ a ) for a < 3m-2 2

Acknowledgments. The first and the third author acknowledge the Post-Doctoral support of École Polytechnique. The third author was also supported by MN grant IdPlus 2011/000661 and by the fellowship START of the Foundation for Polish Science.


Anne Cuzol

Jean-Louis Marchand

Etienne Mémin

Image data assimilation with filtering methods

In this paper we describe several techniques formulated within the stochastic filtering framework for image data assimilation issues. We advocate here the use of hybrid methods between ensemble Kalman methods and particle filters. The former family, despite being theoretically deficient in the sense that it does not in general converge towards the sought-after filtering moments, has demonstrated to be very efficient in practice for high dimensional space filtering issues. At the opposite, the latter are theoretically well posed but face strong practical difficulties in high dimensional spaces. We list here briefly the principal ideas of the underlying hybrid filters, their qualities and their drawbacks. Some comparison results between those different techniques are provided for the filtering of a 2D turbulent flow.

Résumé : Dans cet article nous décrivons plusieurs technique d'assimilation de données images formulées dans le cadre d'un problème de filtrage stochastique non linéaire. Nous prônons l'utilisation de filtres hybrides couplant des filtres de Kalman d'ensemble et les filtres particulaires. La première famille de filtres, bien que déficiente d'un point de vue théorique puiqu'elle ne converge pas vers les moments de la distribution de filtrage cible, a montré son efficacité pour des problèmes d'assimilation de données en très grande dimension. La seconde en revanche, bien posée théoriquement, est confrontée à d'importantes difficultés pratiques en grande dimension. Nous listons brièvement les principes gouvernant la construction de ces filtres, ainsi que leur avantages et défauts. Quelques résultats comparatifs entre ces différentes techniques sont donnés dans le cas du filtrage d'un écoulement turbulent 2D.

Introduction

Data assimilation techniques aim at coupling a system state dynamics with partially observed measurements of this system. Such a procedure is essential for instance in forecasting applications to estimate from observed data a good initial state of the system's variables or to calibrate dynamical parameters. A recent new challenge is emerging in geophysical flow analysis with the availability of high resolution satellite image data [START_REF] Beyou | Weighted ensemble transform kalman filter for image assimilation[END_REF], Corpetti et al., 2009, Titaud et al., 2010]. The assimilation of such data into oceanic or atmospheric numerical models is unfortunately not as easy as it may appear at a very first glance. As a matter of fact, it is immediate to feel like incorporating in weather numerical models complex non linear phenomena such as front, or tornado that are evidently observed in image data sequences, and which are, in the same time, difficult to generate spontaneously from models. The essential difficulty of image assimilation ensues principally from the fact that the observed luminance does not belong to the system state space. In general, there is indeed a complex relation between the luminance function and the flow state variables. This complex relation leads to difficult inverse problems to recover for instance velocity, temperature or pressure variables from a set of images. The other source of concern lies in the fact that images and models do not live at the same scales. As a matter of fact, the spatial resolution of images are nowadays defined on finer spatial grids than models. This still increasing scale difference, due to satellite sensors technological progress, induces a scale relation which complicates further the connection between the state variables and the observations. Besides, large scale models are only imperfectly known. Small scale forcing and geophysical turbulence models for instance must be set up. The alternative, to specify such imperfect large scale evolution models consists in the establishment of stochastic dynamics where the uncertainty is formalized through stochastic processes [START_REF] Franzke | Low-order stochastic mode reduction for a prototype atmospheric gcm[END_REF], Majda et al., 1999, Mémin, 2013]. In this latter case, it is important to point out that we are in a situation where the model random uncertainty is usually much more pregnant than the measurement noise defined at a higher resolution. This extreme case of quasi negligible data noise leads in particular to a badly conditioned situation for data assimilation procedures.

The coupling of a a partially known stochastic dynamics and noisy observations can be expressed through a nonlinear stochastic filtering problem. In the case of image data this filtering problem involves an evolution law of a state variable x(x,t) :

Ω × R + → R d of the form dx(x,t) = M(x(x,t))dt + σ (x,t)dB(t), (1)
accompanied by partial observations of the system through measurements y(x,t) :

Ω × R + → R p y(x,t k ) = H(x(x,t k ) + γ(x,t k ).
(2)

In this system, both the dynamical operator, M, and the observation operator, H, between the measurement space and the state space are, in the general case, nonlinear. The observation equation ( 2) is in addition discrete in time due to a possibly long time interval between two consecutive measurements. The dynamics is at the opposite continuous with a discrete scheme associated to a much shorter time step. Concerning this stochastic dynamics, we will assume the dynamical noise is a Wiener process where the diffusion tensor, σ , may eventually depend on the stochastic process, which gives rise in that case to multiplicative noises. As for the observation noise, γ t k , it can be non-Gaussian, but it is, nevertheless, conditionally independent of the state variable.

Beyond the main difficulty faced here, which concerns essentially the high dimensionality of the system (∼ 10 6 -10 9 ), let us note the modeling of a proper stochastic representation of a state variable dynamics is in general not a straightforward stage. Several different directions can be followed depending on the application targeted. This goes from the definition of stochastic representations of deterministic pdes to the constitution of "inexact" evolution laws incorporating errors or uncertainties. In this article we will not elaborate on such issues and consider that the dynamics has been well defined and follows general form (1).

In the following we will present different solutions devised so far for the filtering of image data. Those techniques have in common to be all based on Monte Carlo principles. The first family of method ensues from a Monte Carlo implementation of the Kalman filtering equations and incorporates efficient procedures to cope with matrix-vector products in high dimension. The second kind of methods concerns the so-called particle filters. It corresponds to a set of methods that approximate the filtering distribution through a linear combination of Dirac masses centered on samples referred to as particles. In the following section we briefly describe those two types of methods. Advantages and deficiencies of both approaches are listed and we then describe, in a subsequent section, several hybrid techniques that aim at combining both families in order to keep the best of their properties. Some comparative results are provided in the last sections.

Assimilation with filtering methods

The first family of filtering methods is intensively used in geophysics. It relies on a Gaussian assumption and extends the Kalman update mechanisms through a Monte Carlo implementation.

Ensemble Kalman filter

The well-known Kalman filter [Kalman, 1960] allows to solve the filtering problem for linear Gaussian state-space models. At each time step, the Gaussian filtering distribution is characterised by its mean and covariance, which are recursively updated through the prediction and correction steps. In order to tackle the filtering problem for non linear and high-dimensional systems, Ensemble Kalman filters have been developed [Evensen, 1994, Houtekamer and Mitchell, 1998, Bishop et al., 2001, Evensen, 2003, Ott et al., 2004]. These methods consist in a Monte Carlo approximation of the filtering distribution with particles that follow the prediction and correction steps of Kalman equations.

At time t k , an approximation of the filtering distribution is given by :

p(x t k |y t 1 :t k ) = N ∑ i=1 δ x (i) t k (x t k ), (3)
where

x (i) t k i = 1, .
. . , N are the particles (also called "ensemble members"). The prediction step consists in propagating the ensemble of particles x a,(i) t k-1 through the non-linear dynamics (including the model noise simulation) in order to obtain the predicted particles, or forecast ensemble {x f ,(i) t k , i = 1, . . . , N}. The empirical ensemble covariance matrix P f e t k is then deduced from the following expression:

P f e k = 1 N -1 N ∑ i=1 x f ,(i) k -x f k x f ,(i) k -x f k T , (4)
where the empirical mean of the forecast ensemble is defined by

x f k = 1 N N ∑ i=1 x f ,(i) k
. The Kalman correction equation is then computed from this covariance matrix. Since the Ensemble Kalman filter is based on a Gaussian assumption, it is known ([Le Gland et al., 2011]) that this filter does not converge towards the true filtering distribution for non linear systems. However, despite this theoretical drawback, this filter has been applied successfully to high-dimensional systems.

Particle filters

Particle filters [START_REF] Gordon | Novel approach to non-linear/non-gaussian bayesian state estimation[END_REF], Doucet et al., 2000[START_REF] Del Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF] are able to solve exactly the filtering equations, up to the Monte Carlo approximation. In this class of methods, the filtering distribution is approximated by a weighted sum of particles :

p(x t k |y t 1 :t k ) = N ∑ i=1 w (i) t k δ x (i) t k (x t k ).
(5) The particle filter is a sequential importance sampling technique: At each time step, particles {x (i) t k , i = 1, . . . , N} are first sampled from an importance distribution π(x t k |x t 0 :t k-1 , y t 1 :t k ) (this corresponds to the prediction step). Then a new weight w (i) t k is computed for each particle (correction step) using the following correction formulae :

w (i) t k ∝ w (i) k-1 p(y t k |x (i) t k )p(x (i) t k |x (i) t k-1 ) π(x (i) t k |x (i) t 0 :t k-1 , y t 1 :t k ) .
The standard choice when applying the particle filter is to fix the proposal distribution to be the transition law p(x t k |x t k-1 ). In that case, it is easy to sample from the proposal (this requires only to sample from the dynamical model). However, this choice has the major drawback of sampling the particles without taking into account the observation y t k . For high-dimensional systems, it has been proved that this standard particle filter can not be efficient [START_REF] Snyder | Obstacles to high-dimensional particle filtering[END_REF].

Different strategies can be adopted in order to make the particle filter more efficient for highdimensional models. First, one can define a precise dynamical model from the observations, aiming at reducing the dimension of the search space. This approach has been used in [START_REF] Avenel | Stochastic level set dynamics to track closed curves through image data[END_REF] for the tracking of closed curves in images sequences. The dynamical model of the curve is constructed from the image data, and the dynamic uncertainty is reduced to a smalldimension noise (corresponding to normal and tangential perturbations of the curve). A low order dynamical representation based on radial basis functions has been used also successfully to track vorticity fields from images [START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF]Memin, 2009]. In both cases, which are associated to random field of low dimension, it has been shown that the simplest version of the particle filter can be successfully used.

Another way to cope with the curse of dimensionality associated to the usual blind sampling strategy of the bootstrap filter (π(x t k |x t 0 :t k-1 , y t 1 :t k ) ∼ p(x t k |x t k-1 )) consists in constructing importance distributions that integrate observations. This approach is presented briefly in the next section.

Hybrid filtering method

Idea

The importance distribution of the particle filter can be constructed in such a way that it integrates the observation y t k at each time step. In particular, one can use the Ensemble Kalman filter described in Section 2.1. This approach, called "Weighted Ensemble Kalman filter", has been proposed in [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]. This filter combines the good properties of the ensemble Kalman mechanics with the correction scheme of the particle filter. Conceptually, this filter is more adapted to a non Gaussian distribution of the ensemble members, which is a situation that may arise at very short time horizon with nonlinear stochastic dynamics. The approach has been extended in [START_REF] Beyou | Weighted ensemble transform kalman filter for image assimilation[END_REF] for the assimilation of image observations, using a specific square-root filter (the Ensemble Transform Kalman filter) for the proposal step. We present below some results obtained with this method.

Application

State-space model In [START_REF] Beyou | Weighted ensemble transform kalman filter for image assimilation[END_REF], the state-space model is described by a stochastic version of the 2D Navier-Stokes equation for incompressible flows :

dξ t = -∇ξ t • wdt + ν∆ξ t dt + ηdB t ( 6
)
where ν is the viscosity coefficient (inverse of the flow Reynolds number) and dB t is a random forcing term. The variable ξ (x) denote the scalar vorticity at point x = (x, y) T , associated to the 2D velocity w(x) = (w x (x), w y (x))

T through ξ (x) = ∂ w y ∂ x -∂ w x ∂ y .
The vector ξ ∈ R |Ω| describes the vorticity over a spatial domain Ω of size |Ω|, and w ∈ R 2|Ω| is the associated velocity field over the same domain.

This dynamical model is associated to a non linear observation equation based on the luminance (color intensity of images) conservation along fluid trajectories, up to a Gaussian noise :

I(x,t k ) = I(x + d(x),t k+1 ) + γ t k (x) ∀x ∈ Ω I , ( 7
)
where I is the luminance function, Ω I is the spatial image domain, and d(x) = t k+1 -δt t k w t (x t )dt denotes the displacement between time t k and time t k+1 , with x k = x. The perturbation γ t k is a Gaussian noise.

Results

The approach has been validated on synthetic turbulence data and tested on real satellite oceanic data.

The first set of experiments concerns a synthetic image sequence with images of size 256 × 256 showing the transport of a passive scalar by a 2D turbulent flow. This sequence is, to some extent, representative of typical satellite images depicting transport processes by oceanic streams, such as sea surface temperature or sea color images. This scalar exhibits a small diffusion, and does not respect strictly a luminance conservation assumption. The 2D turbulence is also maintained by an unknown forcing, which is thus crudely modeled in our dynamical model by a zero mean Gaussian field.

Results are compared to an other data assimilation method which can be used within this framework: the weighted ensemble Kalman filter (WEnKF, [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF]) which also relies on an ensemble Kalman filtering step, integrated within a particle filtering framework. In [START_REF] Papadakis | Data assimilation with the weighted ensemble Kalman filter[END_REF], the WEnKF was used to assimilate linear observations. In case of image data that is nonlinearly related to the state of the system, these linear observations consist in pseudo-observations, i.e. velocity fields (and associated vorticity maps) computed from a given motion estimation technique. In the experiments below, these pseudo-observations are computed from a stochastic version of the well-known Lucas and Kanade motion estimator [START_REF] Corpetti | Stochastic uncertainty models for the luminance consistency assumption[END_REF] applied on each pair of the image sequence, and are called SLK.

The results are first compared in terms of global RMSE between the ground truth vorticity and the estimated vorticity at each time step, and similarly for the motion fields. The WEnKF and WETKF results have been obtained with N = 700 particles. As may be observed on Figure 1-(a), the improvement obtained with the WEnKF filtering technique over SLK measurements is minor, and even deteriorates the results in terms of velocity fields reconstruction. The proposed WETKF technique, directly based on image luminance observations, leads to much better results for the estimation of vorticity and velocity fields.

In addition to global errors comparisons, the analysis of spectra plotted on Figure 1-(b) allows to observe more precisely the accuracy of the different techniques at different scales. As a matter of fact the RMSE constitutes only a performance measure at large scales. It is indeed interesting to note that all methods give relatively coherent results only from the large scales up to the beginning of the inertial range. On the other hand, the result obtained with the WETKF assimilation scheme is closer to the ground truth over the whole scale range (from the largest scales up to the small dissipative scales). This comes from the fact that this technique can take benefit of all available information in the image data, while pseudo-observations given by a local motion estimator will not be able to improve the estimation at all scales. Filters based on pseudo-observations are unable to correct the loss of energy caused by the smoothing operator used in the external estimation procedure.

For a qualitative visual comparison, estimated vorticity maps are plotted on Figure 2 satellite images of Sea Surface Temperature (SST), with large areas of missing data due to the cloud cover and presence of coastal regions. The sequence consists of 48 images of size 256×256.

The images have a spatial resolution of 0.1 degree (10 km) and a temporal latency of 24 hours. The sequence is centered on an area of the Pacific ocean off the Panama isthmus and shot during an El Niño-Southern oscillation. Representative results at different time instants of the sequence are shown in Figure 3. The first column of Figure 3 shows the estimated velocity fields at times k = 1, 10, 24 and k = 39, superimposed on the corresponding SST images. The second column shows the estimated vorticity maps and velocity fields. The initialization for image 1 is based on the estimation provided by the local motion estimation approach (SLK). We see that this initialization provides only a rough large scale motion field. This estimate is refined afterward by the filtering process. We can note that the motion fields estimated along the sequence stick quite well to the big image structures observed on the SST images. The sequence of motion fields does not seem to be perturbed by the big missing data regions observed intermittently. Finally, we can note that the result obtained shows well an increase of the turbulent agitation with an intensification of the El-Niño phenomenon. However, one can identify many small scale vortices in the estimation results that may not be physically relevant. A better characterization of the dynamical noise could help correcting these artefacts.

Conditional simulation for image data assimilation

The image data (given by satellites for instance) offer a high spatial resolution but a low temporal resolution, so that the time step between observations can be long in the assimilation scheme. This may cause discontinuities in the assimilation result, since the filtering relies on the dynamical model only between observations. In order to reduce such discontinuities and avoid the estimation of state trajectories that are not physically consistent, one can use a Monte Carlo sequential smoothing technique based on conditional simulation of diffusions [START_REF] Delyon | Simulation of conditioned diffusions and applications to parameter estimation[END_REF]. This technique, proposed in [START_REF] Cuzol | Monte carlo fixed-lag smoothing in state-space models[END_REF], relies on the filtering result and allows the sampling of continuous smoothed trajectories between observation times. The method can be applied to non linear and multidimensional models. An example is given in Figure 4, related to the estimation of a turbulent flow from images. The experimental setup is similar to the one presented in Section 2.3. The filtering and smoothing results are given as illustration between two given observation times t k and t k+1 . The temporal discontinuity in the filtering result can be seen on the first row, while the second row shows the result after smoothing. This promising technique can be extended to build alternative data assimilation techniques for data that are observed without noise (or very small noise). In that case, the assimilation problem can be formulated as the conditional simulation of a partially observed model, where the state (the vorticity for instance) is gradually corrected toward the observations (image data). The aim is then to avoid the discontinuities between observation times and to provide smoothed trajectories in a sequential way without post-processing. This work in progress is based on the theoretical results presented in [Marchand, 2012, Marchand, 2013].

Conclusion and perspectives

In this paper we have shown how a Kalman ensemble technique can be embedded within a particle filter. The resulting technique compared to ensemble filtering techniques allows robustifying the filtering and provides improved results. It opens also a way to assimilate different observations attached either to the particle filtering stage or to the sampling step, encoded through the ensemble Kalman filter mechanism. Such an ability has not been yet experimentally assessed and we plan to explore it in a near future. The hybrid filters show standard limitations associated to filtering issues where a continuous dynamics is coupled with a discrete sequence of observations. In that case an implausible trajectory is generated by a violent correction at observation times. We show briefly how to correct a posteriori this deficiency with a fixed lag smoothing. This conditional simulation approach can be extended to built a data assimilation method associated to the ill-posed case of low observation noise. We are currently exploring this promising technique.

FIGURE 1 .

 1 FIGURE 1. Synthetic 2D turbulent flow. (a) RMSE between mean estimate of vorticity and ground truth; (b) Energy spectra.

  at a given time instant. The scalar image observation is first presented, together with the ground truth vorticity on Figure 2(b).

FIGURE 2 .

 2 FIGURE 2. Synthetic 2D turbulent flow. (a) Example of scalar image of the sequence at given time k; (b) Ground truth vorticity at time k; (c) SLK vorticity estimate; (d) Mean estimate of vorticity with WEnKF assimilation of SLK observation; (e) Mean estimate of vorticity with WETKF direct assimilation of image data.

  Soumis au Journal de la Société Française de Statistique File: SFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 14, 2013

  Soumis au Journal de la Société Française de Statistique File: SFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 14, 2013

  Soumis au Journal de la Société Française de Statistique File: SFDS.tex, compiled with jsfds, version : 2009/12/09 date: October 14, 2013