Edwige Godlewski
email: edwige.godlewski@upmc.fr

Martin Parisot
email: martin.parisot@inria.fr

Jacques Sainte-Marie
email: jacques.sainte-marie@inria.fr

Fabien Wahl
email: fabien.wahl@upmc.fr

CONGESTED SHALLOW WATER MODEL: ROOF MODELING IN FREE SURFACE FLOW

Keywords: May 2, 2018. 2010 Mathematics Subject Classification. 35L65, 76B07, 76M12, 35Q35, 74F10 Shallow water equations, Congested hyperbolic model, Unilateral constraint, Well-balanced scheme, Entropic scheme

We are interested in the modeling and the numerical approximation of flows in the presence of a roof, for example flows in sewers or under an ice floe. A shallow water model with a supplementary congestion constraint describing the roof is derived from the Navier-Stokes equations. The congestion constraint is a challenging problem for the numerical resolution of hyperbolic equations. To overcome this difficulty, we follow a pseudo-compressibility relaxation approach. Eventually, a numerical scheme based on a Finite Volume method is proposed. The well-balanced property and the dissipation of the mechanical energy, acting as a mathematical entropy, are ensured under a non-restrictive condition on the time step in spite of the large celerity of the potential waves in the congested areas. Simulations in one dimension for transcritical steady flow are carried out and numerical solutions are compared to several analytical (stationary and non-stationary) solutions for validation.

Introduction

We deal with the derivation and the numerical resolution of a model describing shallow water flows in the presence of a roof. The roof denotes an impermeable surface above the flow which constrains the water surface in several parts of the domain. We do not assume that the roof is fixed, but the roof has a prescribed motion and the surface pressure is seen as an unilateral constraint. This assumption is not really essential to the analysis but greatly simplifies the reading. In particular the numerical procedure is not affected by the roof dynamics. The description of the fluid/structure interaction is studied in the literature for the Navier-Stokes equations [START_REF] Quarteroni | Computational vascular fluid dynamics: problems, models and methods[END_REF]. However, for vertical-integrated models such as the shallow water model, the interaction with a roof has been considered only recently. The main applications of these models concern the prediction of urban flooding taking into account the flows in sewers or the hydraulic network in small hydroelectric power plants, [START_REF] Capart | Numerical and experimental water transients in sewer pipes[END_REF][START_REF] Fuamba | Contribution on transient flow modelling in storm sewers[END_REF]. Moreover, the production of sustainable energy from wave power using floating buoys can also be realized considering a dynamical roof.

For the modeling of partially free surface flow, i.e. flow where some parts of the domain have a free surface and do not reach the roof and some other parts are in contact with the roof, two main approaches are proposed in the literature.

The first one can be classified as a coupling strategy between the free surface and the pressurized areas [START_REF] Kashiwagi | Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method[END_REF][START_REF] Lannes | On the dynamics of floating structures[END_REF]. The main advantage of this approach relies on the computation of the surface pressure induced by the roof using an elliptic equation obtained thanks to Bernoulli's equation. In addition, each sub-domain can be solved with the most appropriate numerical strategy. However, this strategy assumes a smooth enough flow in order to write transmission conditions at the coupling interface. This assumption is not always satisfied. In particular, it is wellknown that hydraulic models can lead to discontinuous free surfaces which imply an advection of a discontinuous pressure. In addition, this strategy requires the description of the interface position which is a challenging issue since the interface is strongly linked to the flow dynamics.

The second strategy to model partially free surface flow is based on the design of a unified model in the whole domain. Our work is based on this idea. Previous works also deal with a unified model [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF][START_REF] Fuamba | Contribution on transient flow modelling in storm sewers[END_REF]. This strategy is advantageous because it does not require the description of the interface position and the transmission conditions at the interface. The partially free surface models proposed in the cited works can be described as shallow water models with an additional pressure given as an explicit function of the water depth. The additional pressure comes from the compressibility of the fluid in the pressurized areas and introduces the sound celerity in the water. The material velocity is very small in comparison to the sound celerity which leads to the particular low-Mach regime. However, in our approach, we do not consider the water as compressible. In Section 1.1, a shallow water model with a constraint, from now on called congestion constraint, modeling the roof is derived from the Navier-Stokes equations in a tank under the classical shallow water assumption [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF].

A consequence of this approach is that the model we obtain is of hyperbolic type with a supplementary congestion constraint. This type of equations appears in many contexts. A congested hyperbolic model describing pressureless gas dynamics is proposed in [START_REF] Bouchut | A hierarchy of models for twophase flows[END_REF]. The model is similar to a model of crowd motion [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF]. In [START_REF] Berthelin | A traffic-flow model with constraints for the modeling of traffic jams[END_REF], the authors use a congested hyperbolic model to describe traffic jams. In [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF], the authors use a degenerated parabolic model to describe porous media. For some particular parameters, the model can be written as a congested hyperbolic model. In [START_REF] Perrin | Free/congested two-phase model from weak solutions to multi-dimensional compressible Navier-Stokes equations[END_REF], a congested hyperbolic model is presented as the singular limit of a barotropic model for granular flow. One can also highlight the relevance of congestion constraints in models for chemotaxis in order to fix mathematical blowup [START_REF] Perthame | PDE models for chemotactic movements: Parabolic, hyperbolic and kinetic[END_REF].

One of the first analysis of a congested problem is realized in [START_REF] Lions | On a free boundary barotropic model[END_REF]. In [START_REF] Berthelin | Weak solutions for a hyperbolic system with unilateral constraint and mass loss[END_REF] a congested hyperbolic system where the constraint is verified by mass loss is analyzed. In a more general framework, congested linear hyperbolic conservation laws are analyzed in [START_REF] Després | A geometrical approach to nonconservative shocks and elastoplastic shocks[END_REF][START_REF] Després | Weak solutions to Friedrichs systems with convex constraints[END_REF]. Quasi-linear scalar conservation laws are studied in [START_REF] Levi | Obstacle problems for scalar conservation laws[END_REF]. The congested pressureless model introduced in [START_REF] Bouchut | A hierarchy of models for twophase flows[END_REF] is analyzed in [START_REF] Berthelin | Existence and weak stability for a pressureless model with unilateral constraint[END_REF]. Compared to the previous models, the model introduced in our work presents a bound of congestion depending on space and time.

The numerical resolution of the above cited models [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF][START_REF] Fuamba | Contribution on transient flow modelling in storm sewers[END_REF] is a challenging problem due to the congestion constraint. In particular in our context, where the upper bound of the water level, called the opening in the sequel, is a function depending on time and space. In Section 1.2, a strictly hyperbolic approximate model is proposed which is based on a pseudo-compressibility approximation. The approximate model is similar to the above cited models [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF][START_REF] Fuamba | Contribution on transient flow modelling in storm sewers[END_REF] with the difference that the relaxation parameter is not linked to any physical process. Nevertheless, we will require a similar asymptotic limit, namely the high potential force regime [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF].

In addition, since the opening is a space (and time) function, a particular attention is given to the well-balanced property, i.e. the numerical preservation of the lake at rest. We emphasize that this property is also required in the congested area to get a relevant approximation of the surface pressure.

The main achievement is the description of a relaxation strategy, see Section 1.2 and the adaptation of a numerical scheme for the approximation of the congested shallow water model, see Section 2. Indeed, the numerical scheme for an advection model in high potential force regime [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF] can be easily adapted to our (or other) congested hyperbolic models. The well-balanced property and the dissipation of the mechanical energy, acting as a mathematical entropy, are proven. Eventually, several numerical simulations in a one dimensional framework are performed and compared to analytical solutions. Simulations of discontinuous solutions are also proposed.

Mathematical modeling

1.1. Formal derivation of the congested shallow water model. A shallow water type model with an additional congestion constraint modeling a roof, i.e. an impermeable boundary above the flow, is derived from the Navier-Stokes equations using classical thin layer arguments. Let us consider a flow contained between two given surfaces respectively called bottom and roof. We assume that the two surfaces can be parametrized by two mono-valued smooth enough functions B (x, t) (for the bottom) and R (x, t) (for the roof) which satisfy B (x, t) ≤ R (x, t). We consider the domain Ω The bottom (resp. the roof) moves according to its own prescribed velocity with

(t) = Ω x × [B (x, t) , R (x, t)] where Ω x ⊂ R d-1 with d = 2 or d = 3, see Figure 1. The space variable is split into its horizontal component x ∈ Ω x and its vertical component z ∈ [B (x, t) , R (x, t)].
U B (x, t) ∈ R d-1 (resp. U R ) the horizontal component and W B (x, t) ∈ R (resp.
W R ) the vertical component of the given velocity. More precisely for Γ = B or Γ = R, it satisfies

∂ t Γ + U Γ • ∇ x Γ -W Γ = 0. (1.1)
The fluid is assumed homogeneous, with a density set to 1 for readability. The part of the domain occupied by the fluid, denoted Ω f (t), is marked by the color function ψ (x, z, t) ∈ {0, 1} such that ψ (x, z, t) = 1 iff (x, z) ∈ Ω f (t). More precisely we assume that the fluid surface is regular enough and that there exists a mono-valued function η (x, t) parametrizing the free surface such that ψ writes

ψ(x, z, t) = 1, if B(x, t) ≤ z ≤ η(x, t), 0, otherwise.
The velocity of the fluid satisfies the incompressible Navier-Stokes equations with gravity, i.e. for any (x, z, t)

∈ Ω f (t) × R +    ∇ x • u + ∂ z w = 0, ∂ t u + (u • ∇ x ) u + w∂ z u + ∇ x p -∇ x • σ xx -∂ z σ xz = 0, ∂ t w + u • ∇ x w + w∂ z w + ∂ z (p + gz) -∇ x • σ zx -∂ z σ zz = 0, (1.2)
with u (x, z, t) ∈ R d-1 the horizontal component and w (x, z, t) ∈ R the vertical component of the velocity, p (x, z, t) ∈ R the pressure, g the gravity acceleration and σ (u, w) the viscosity tensor. We denote by n η the outward normal to the fluid surface defined by

n η = 1 1 + |∇ x η| 2 -∇ x η 1 .
The fluid is advected by the flow and the color function describing the part of the domain occupied by the fluid satisfies the following relation

∂ t ψ + u • ∇ x ψ + w∂ z ψ = 0. (1.3)
From now on, a Newtonian fluid is considered, i.e. the viscosity tensor is given by

σ = σ xx σ xz σ zx σ zz with σ xx = 2µD x (u) , σ xz = µ (∂ z u + ∇ x w) , σ zz = 2µ∂ z w, σ zx = µ (∂ z u + ∇ x w) t ,
where

D x (u) = 1 2 ∇ x u + (∇ x u)
t is the symmetric gradient and µ > 0 the viscosity coefficient. For each surface Γ = B or Γ = R the outward unit normal, see Figure 1, denoted by n Γ , is defined by

n B = 1 1 + |∇ x B| 2 -∇ x B 1 and n R = 1 1 + |∇ x R| 2 ∇ x R -1 .
At the bottom and at the roof, a no-penetration condition is considered, i.e.

ψ |z=R u |z=R -U R w |z=R -W R • n R ≥ 0 and u |z=B -U B w |z=B -W B • n B = 0. (1.4)
The inequality at the roof is required to allow the fluid to detach itself from the surface. In addition, the friction between the fluid and the surface Γ = B or Γ = R is modeled by the following Navier condition

ψ |z=Γ Π Γ (σ -pI d ) n Γ -κ Γ u |z=Γ -U Γ w |z=Γ -W Γ = 0 (1.5)
with κ Γ > 0 a friction coefficient at the surface Γ and Π Γ the projection on the tangential plane to the surface. At the free surface, the continuity of the stresses reads

1 -ψ |z=R (σ + (P -p) I d ) n η = 0 (1.6)
with P (x, z, t) a given atmospheric pressure. At this stage of modeling, when the pressure in the fluid reaches the atmospheric pressure, i.e. p = P , the fluid can either create vacuum, i.e. ψ = 0, or the pressure can decrease below the atmospheric pressure p < P , see [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF] for more details. In [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF] the authors propose a numerical strategy able to forbid the creation of vacuum but it is not clearly described at the continuous level. In the current work, only non-negative relative pressure is considered, i.e.

p |z=η ≥ P |z=η .
(1.7)

The drawback of this assumption is that areas without water trapped under the roof can appear, which is not physically relevant. Some perspectives to overcome this drawback are discussed in the conclusion section.

The system (1.1)-(1.7) has to be completed with compatible initial data. More precisely the fluid domain at the initial time reads Ω f (0) = Ω 0 f with Ω 0 f given and η (x, z, 0) = η 0 (x). Similarly, for any (x, z) ∈ Ω 0 f , we assume u (x, z, 0) = u 0 (x, z) and w (x, z, 0) = w 0 (x, z) where u 0 and w 0 are given. Note that to be physically relevant, the initial fluid has to be contained in the domain, i.e. Ω 0 f ⊂ Ω (0) and the velocity should satisfy the incompressibility condition, i.e. first equation of (1.2).

The following domain invariance holds.

Proposition 1.1. If Ω 0 f ⊂ Ω (0), then Ω f (t) ⊂ Ω (t) for any time t.

Proof. By integration of (1.3) in R d \ Ω (t) and using a Leibniz integral rule, we get

∂ t R d \Ω(t) ψ dx dz = R d-1 ψ |z=B u |z=B -U B • ∇ x B -w |z=B -W B dx - R d-1 ψ |z=R u |z=R -U R • ∇ x R -w |z=R -W R dx.
Using the non-penetration inequality (1.4) and integrating in time, it yields

R d \Ω(t) ψ (x, z, t) dx dz ≤ R d \Ω(0) ψ (x, z, 0) dx dz.
The assumption on the initial condition allows to conclude since ψ (x,

z, t) = 1 iff (x, z) ∈ Ω f (t).
This property ensures that the fluid stays confined between the bottom and the roof.

1.1.2. Congested shallow water model. Let us introduce the congested shallow water model we are interested in. The unknowns of the model are the water depth h, the vertical-averaged horizontal velocity u and the surface pressure p η which satisfy

     ∂ t h + ∇ x • (hu) = 0 , ∂ t (hu) + ∇ x • hu ⊗ u + g 2 h 2 I d = -h∇ x (gB + p η ) -κ (u -U ) +2µ (∇ x • (hD x (u)) + ∇ x (h∇ x • u)) , (1.8)
with the following given parameters: κ the efficient friction parameter, U the efficient velocity of forcing and µ the viscosity coefficient. In addition, the water depth has to satisfy the congestion constraint

h ≤ H (1.9)
with the opening between the roof and the bottom H (see Figure 1) defined by

H (x, t) = R (x, t) -B (x, t) .
Eventually, the surface pressure has to satisfy [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF]. More precisely, considering the two characteristic space dimensions along the vertical and the horizontal dimensions, respectively H and L, and the characteristic time of observation T , the dimensionless variables read û =

T L u, ŵ = T H w, p = T L 2 p, ĝ = H U 2 g, κΓ = T L κ Γ and μ = T L 2 µ.
For readability, the hat ˆis dropped even if from now on the variables are dimensionless. We claim the following proposition. Proposition 1.2. For ε = H L small enough, assume that the parameters satisfy the following scaling

∇ x B = O (ε) , ∇ x R = O (ε) , ∇ x η = O (ε) µ = εµ 0 , κ B = εκ 0
B , and κ R = εκ 0 R . Then (1.8)-(1.10) with the initial conditions

h 0 (x) = η 0 (x) -B (x, 0) and u 0 (x) = 1 h 0 (x) η 0 (x) B(x,0) u 0 (x, z) dz
and the parameters given by

κ = κ B + 1 h≥H κ R , U = 1 κ κ B U B + 1 h≥H κ R U R , µ = εµ 0 and P η (x, t) = P (x, η (x) , t) (1.11) with κ R =              κ 0 R 1 + εh 3µ 0 1 h≥H κ 0 R -1 2 κ 0 B , if κ 0 B ≤ 2κ 0 R κ 0 R 1 - εh 3µ 0 1 h≥H κ 0 R - 1 2 κ 0 B , otherwise and κ B =              κ 0 B 1 + εh 3µ 0 κ 0 B -1 h≥H 1 2 κ 0 R , if κ 0 R ≤ 2κ 0 B κ 0 B 1 - εh 3µ 0 κ 0 B -1 h≥H 1 2 κ 0 R , otherwise
is formally an approximation in O ε 2 of (1.1)-(1.7). More precisely

h (x, t) -(η (x, t) -B (x, t)) = O ε 2 and u (x, t) - 1 h (x, t) η(x,t) B(x,t) u (x, z, t) dz = O ε 2 .
Proof. As it was done in [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF], the vertical integration of (1.2) yields the mass conservation (first equation of (1.8)) and the following momentum balance

∂ t η B u dz + ∇ x • η B (u ⊗ u + pI d ) dz = p |z=η ∇ x η -p |z=B ∇ x B -κ B u |z=B -U B -1 h≥H κ R u |z=R -U R + ∇ x • η B σ xx dz .
(1.12) Thanks to Proposition 1.1, ψ |z=R (x, t) = 1 h≥H (x, t). Let us start by formally identifying a first order approximation. As in [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF], the main term of the third equation of (1.2) reads p = p |z=η + g (η -z) + O (ε). The main term of the second equation of (1.2) and the Navier boundary conditions (1.5) lead to u = u + O (ε). Injected in (1.12), the momentum balance of the congested shallow water model reads at first order

∂ t (hu) + ∇ x • hu ⊗ u + g 2 h 2 I d = -h∇ x gB + p |z=η -κ 0 B (u -U B ) -1 h≥H κ 0 R (u -U R ) .
In order to get a better approximation, the first perturbation terms are taken into account. The vertical momentum equation of (1.2) leads to

p = p η + g (η -z) -2εµ 0 ∇ x • u + O ε 2
where the surface pressure reads p η = p |z=η -2εµ 0 ∂ z w |z=η . Considering the horizontal momentum equation of (1.2) and the Navier conditions at the bottom and the roof (1.5), the horizontal velocity of (1.2) satisfies

             ∂ 2 zz u = - ε µ 0 h κ 0 B (u -U B ) + 1 h≥H κ 0 R (u -U R ) + O ε 2 , ∂ z u |z=B = εκ 0 B µ 0 (u -U B ) + O ε 2 , 1 h≥H ∂ z u |z=R = -1 h≥H εκ 0 R µ 0 (u -U R ) + O ε 2 .
By taking into account the definition of u, one can check that the horizontal velocity of (1.2) is under the form

u = α 1 2 (z -B) 2 + α 2 (z -B) + α 3 + O ε 2 with α 1 = - ε µ 0 h κ 0 B (u -U B ) + 1 h≥H κ 0 R (u -U R ) , α 2 = ε µ 0 κ 0 B (u -U B ) ,
and

α 3 = u + εh µ 0 1 h≥H κ 0 R 6 (u -U R ) - κ 0 B 3 (u -U B ) .
Inserted in (1.12), the shallow water model (1.8) follows, provided we define U and κ by (1.11).

It remains to establish (1.9) and (1.10). A direct consequence of Proposition 1.1 leads to B ≤ η ≤ R and thus to the congestion constraint (1.9). Similarly, the non-negativity of the relative pressure (1.7) implies p η ≥ P |z=η . Introducing the first equation of (1.6) into the second equation of (1.6) gives

(1 -1 η≥R ) p η -P |z=η = O ε 2 .
We get the expected constraint (1.10).

The two cases in the definition of κ B and κ R are chosen to get positive parameters whatever the water depth. More precisely, both correspond to the same approximation at order O ε 2 . It follows that the effective friction parameter κ is always positive.

1.1.3. Preliminary analysis of the congested shallow water model. Let us mention some mathematical properties of the congested shallow water model. The objective is to highlight the physical relevance of the model (1.8)-(1.10). For a deeper analysis, we refer to [START_REF] Lannes | On the dynamics of floating structures[END_REF]. In particular, the dissipation of mechanical energy is a crucial point for many applications such as renewable energy production using buoys. From the mathematical point of view, the energy balance, acting as a mathematical entropy, is an argument for the existence of long time solutions. In addition the conservation of the steady state, in particular the lake at rest, is relevant for many applications because the flow is close to this state. To characterize these two properties, the potential of the conservative forces defined by

φ = g (h + B) + P η (1.13)
is a useful quantity.

∂ t E + ∇ x • G = -(p η -P η ) ∂ t H + h∂ t (gB + P η ) - κ 2 u 2 -U 2 + u -U 2 -2µh D x (u) 2 + |∇ x • u| 2
with the mechanical energy E = K + P, the kinetic energy and the potential energy respectively given by

K = 1 2 h u 2 and P = gh B + h 2 + hP η
and the energy flux

G = 1 2 u 2 + g (h + B) + p η -2µ (∇ x • u + D x (u)) hu.
Proof. The kinetic energy balance law is obtained by multiplying the momentum balance of (1.8) by the mean velocity u which gives

∂ t K + ∇ x • (K -2µh (∇ • u + D x (u)) u) = -hu • ∇ x (φ + p η -P η ) - κ 2 u 2 -U 2 + u -U 2 -2µh D x (u) 2 + |∇ x • u| 2
with φ defined by (1.13). Furthermore, the potential energy balance law is obtained by multiplying the mass conservation of (1.8) by the potential φ. It yields

∂ t P + ∇ x • ((φ + p η -P η ) hu) = hu • ∇ x φ + ∇ x • ((p η -P η ) hu) + h∂ t (gB + P η ) .
Summing the two energy balance laws, it remains to estimate the source term

(p η -P η ) ∇ x • (hu) = -(p η -P η ) ∂ t h -H -(p η -P η ) ∂ t H.
Eventually, considering the constraint (1.10), we get

0 = ∂ t (p η -P η ) h -H = (p η -P η ) ∂ t h -H + h -H ∂ t (p η -P η ) which implies that (p η -P η ) ∂ t h -H = 0 since p η -P η = 0 or h -H = 0.
The congested shallow water model (1.8)-(1.10) admits the following lake at rest solution, i.e. a steady solution with a vanishing velocity. To provide an expression of it, we use the positive part function defined by (ψ

) + = |ψ| + ψ 2
and we consider a uniform potential of conservative forces φ 0 .

Proposition 1.4. Let ∂ t R = ∂ t B = ∂ t P η = U = 0.
Then for any φ 0 ∈ R, the lake at rest defined by

u = 0 and h = min φ0-Pη g -B + , H (1.14)
is a steady solution of (1.8)-(1.10) and the surface pressure reads

p η = P η + 1 φ0≥g(B+H)+Pη (φ 0 -gR -P η ) .
Proof. The lake at rest solution defined by (1.14) satisfies

h (g (B + h) + p η -φ 0 ) = 0.
In addition, since u = 0, we conclude that ∂ t (hu) = 0 and ∂ t h = 0. Eventually, it is clear that the lake at rest (1.14) ensures the congestion constraint (1.9) and (1.10).

The lake at rest is not the only steady solution of (1.8)-(1.10). However, it is the simplest one. Some works in the literature consider more complex steady states [START_REF] Castro | Well-balanced high order extensions of Godunov's method for semilinear balance laws[END_REF][START_REF] Jin | Two interface-type numerical methods for computing hyperbolic systems with geometrical source terms having concentrations[END_REF][START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography or manning friction[END_REF][START_REF] Noelle | High-order well-balanced finite volume WENO schemes for shallow water equation with moving water[END_REF].

1.2. Relaxation of the congestion constraint. Due to the congestion constraint (1.9), the congested shallow water model is difficult to handle analytically as well as numerically. For this reason, we propose new formulations of the congestion constraint. In the context of incompressible fluid dynamics, several strategies were proposed in the literature to deal with a divergence free constraint. These strategies require that the constraint is written under the form of a divergence. Thus let us first write the unilateral constraint of the congested shallow water model under the form of a divergence. In such a context, the surface pressure p η acts as a Lagrange multiplier in the congested areas.

Proposition 1.5. Assume that h is smooth enough and h 0 (x) ≤ H (x, 0).

Then the following constraints

i) h ≤ H ii) 1 h≥H ∂ t H + ∇ x • (hu) = 0 iii) ∂ t h -H + = 0
are equivalent for any solution of (1.8) and (1.10).

Proof. Let (h, u, p η ) be solution of (1.8) and (1.10).

• i)⇒ii): Assume i). Multiplying the mass conservation equation of (1.8) by the indicator function for h ≥ H, the constraint ii) is recovered since in the congested area we have h = H. • ii)⇒iii): Assume ii). Using the mass conservation, the constraint becomes

1 h≥H ∂ t h -H = 0. • iii)⇒i): Assume iii). Since h 0 (x) ≤ H (x, 0), iii) leads to i). Remark 1.
Once the congestion constraint is under the divergence form ii) of Proposition 1.5, one can remark that in the congested area, i.e. h = H, the surface pressure satisfies

-∇ x • (h∇ x p η ) = ∇ x • ∇ x • hu ⊗ u + g 2 h 2 I d + gh∇ x B +∇ x • (κ (u -U )) -∂ 2 tt H -2µ∇ x • (∇ x • (hD x (u)) + ∇ x (h∇ x • u)) .
(1.15)

In order to recover the surface pressure from (1.15), boundary conditions are required. In [START_REF] Lannes | On the dynamics of floating structures[END_REF], the continuity of the surface pressure is assumed. However, this strategy assumes a smooth enough flow, which is not always satisfied since hyperbolic models can lead to discontinuous solutions, see Section 2.2.3. In addition, this strategy requires the description of the interface position which is a challenging issue since the interface is strongly linked to the flow dynamics.

The relaxation method is a strategy frequently used to approximate the divergence free constraint, see [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF][START_REF] Témam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (I)[END_REF]. We introduce λ > 0 the relaxation parameter, characterizing the amplitude of the penetration in the roof. Let us write the approximate model

     ∂ t h λ + ∇ x • (h λ u λ ) = 0, ∂ t (h λ u λ ) + ∇ x • h λ u λ ⊗ u λ + g 2 h 2 λ I d = -h λ ∇ x (gB + p λ ) -κ (u λ -U ) +2µ (∇ x • (hD x (u λ )) + ∇ x (h∇ x • u λ )) (1.16) completed with the initial conditions h λ (x, 0) = h 0 λ (x) and u λ (x, 0) = u 0 λ (x)
. The unknowns of (1.16) are h λ the approximated water depth which does not exactly satisfy (1.9) and u λ the approximated vertical-averaged horizontal velocity. The unknown surface pressure p η in (1.8) is replaced by p λ which is taken as an implicit function of h λ . More precisely, we assume

λ 2 L (p λ -P η ) = g1 h≥H ∂ t H + ∇ x • (hu) .
(1.17)

We assume that the relaxation operator L : Λ → L (Λ) is invertible. However since choosing, for example, L = -∆ x needs to specify a boundary value problem, we will not discuss the domain of definition and the domain of invertibility of the operator L in the general case and we start by discussing the choice of L. Classically the regularization operator L is chosen among [START_REF] Shen | Pseudo-compressibility methods for the unsteady incompressible Navier-Stokes equations[END_REF] where these operators are presented in the context of Navier-Stokes equations. However, this list is not exhaustive. Since the congestion constraint can be written under the form of a time derivative, see iii) of Proposition 1.5, it seems interesting to choose L under the form

I d , -∆ x , ∂ t or -∂ t ∆ x , see
L (Λ) = ∂ t L (Λ)
with L an operator chosen among I d or -∆ x . Then L is invertible and L -1 is a strictly positive map. We will however skip the theoretical details and do not precise the spaces on which this property holds when L = I d . Note that with this choice, an additional initial condition on the surface pressure p λ (x, 0) = p 0 λ (x) is required and is not provided by (1.8). However, the equation (1.17) can be integrated in time. More precisely choosing the initial condition on p λ such that

λ 2 L p 0 λ (x) -P η (x, 0) = g h 0 λ (x) -H (x, 0) + the surface pressure reads for any (x, t) ∈ R d-1 × R + λ 2 L (p λ (x, t) -P η (x, t)) = g h λ (x, t) -H (x, t) + .
(1.18)

This choice corresponds to a pseudo-compressible approximation, see [START_REF] Shen | Pseudo-compressibility methods for the unsteady incompressible Navier-Stokes equations[END_REF]. In the context of incompressible fluid dynamics, the convergence of the approximate model to the constrained model when λ goes to 0 using the classical regularization operators cited above, can be mathematically proven and the error introduced by the relaxation can be estimated, see for example [START_REF] Rannacher | On Chorin's projection method for the incompressible Navier-Stokes equations[END_REF].

The regularization operator L should be chosen in order to preserve the lake at rest solution, see Proposition 1.4 as well as the positivity of the relative pressure (1.10). In addition, the regularization operator should ensure the dissipation of the mechanical energy given in Proposition 1.3. In the following the latter properties are proven in the case of pseudo-compressible approximation (1.18). To characterize these properties, the potential of the conservative forces which now writes (see (1.13))

φ λ = g (h λ + B) + p λ is again introduced. Let us start by a result on the energy balance. We skip the proof which is similar to the proof of Proposition 1.3. Proposition 1.6. For smooth enough solutions of (1. [START_REF] Chanson | Hydraulics of Open Channel Flow[END_REF]) and (1.18), the following dissipation law holds

∂ t E λ + ∇ x • G λ = -(p λ -P η ) ∂ t H + h λ ∂ t (gB + P η ) - κ 2 u λ 2 -U 2 + u λ -U 2 -2µh λ D x (u λ ) 2 + |∇ x • u λ | 2
with the mechanical energy E λ = K λ + P λ , the kinetic energy and the potential energy respectively given by

K λ = 1 2 h λ u λ 2 and P λ = gh λ B + h λ 2 + h λ P η + R (h λ ) ,
the potential energy relative to the regularization R (ψ) such that

d ψ R = L -1 g λ 2 ψ -H + and the flux of energy G λ = 1 2 u λ 2 + g (h λ + B) + p λ -2µ (∇ x • u λ + D x (u λ )) h λ u λ .
Let us now prove the positivity of the relative pressure (1.10).

Proposition 1.7. The solution of (1. [START_REF] Chanson | Hydraulics of Open Channel Flow[END_REF]) and (1.18) satisfies

1 -1 h λ ≥H (p λ -P η ) = 0 and p λ ≥ P η .
Proof. The positivity of the relative pressure is easy to obtain from the expression (1.18) of the surface pressure. Now multiplying (1.18) by the indicator function of the areas where the water does not reach the roof yields

1 -1 h λ ≥H (p λ -P η ) = 1 -1 h λ ≥H L -1 g λ 2 h λ -H + .
We conclude since L -1 is a strictly positive map.

Eventually, the conservation of the lake at rest is highlighted.

Proposition 1.8. Assume ∂ t R = ∂ t B = ∂ t P = U = 0.
Then for any φ 0 ∈ R, the lake at rest defined by u λ = 0 and

h λ =        φ 0 -P η g -B + , if φ 0 ≤ g H + B + P η , I d + λ 2 L -1 H + λ 2 L φ 0 -P η g -B , else (1.19)
is a steady solution of (1.16) and (1.18) and

p λ = max (φ 0 -g (B + h λ ) , P η ) . (1.20)
Proof. Let us first consider the case where φ 0 > g H + B + P η . Starting from h λ + λ 2 L (h λ ) using (1.19), we get

λ 2 L (φ 0 -g (B + h λ ) -P η ) = g h λ -H .
Using (1.18) it leads to

λ 2 L (φ 0 -g (B + h λ ) -P η ) = λ 2 L (p λ (x, t) -P η (x, t))
and composing by L -1 we get p λ = φ 0 -g (B + h λ ) then φ λ = φ 0 .

In the case where φ 0 ≤ g H + B + P η , (1.19) leads to φ 0 = g (h λ + B) + P η . Moreover (1.18) leads to p λ = P η . Thus we finally get φ λ = φ 0 .

In addition, since u λ = 0, we conclude that ∂ t (h λ u λ ) = 0 and ∂ t h λ = 0.

At the limit λ → 0, one can check that the state (1.19) tends to the state (1.14). It is a (very weak) argument to claim that (1.16) and (1.18) is a consistent approximation of (1.8)-(1.10). It is consistent (as λ → 0) at least for some steady solutions.

To conclude, the pseudo-compressibility method seems a suitable approximation of the congested shallow water model (1.8)-(1.10). In particular, the simplest pseudo-compressibility method leading to an approximation satisfying the properties of the congested shallow water model (1.8)-(1.10) is given by L = I d . In addition, in this case, (1.18) gives an explicit function for the surface pressure

p λ (h λ ) = P η + g h λ -H + λ 2 . (1.21) Substituting it in (1.16) leads to                  ∂ t h λ + ∇ x • (h λ u λ ) = 0 , ∂ t (h λ u λ ) + ∇ x • h λ u λ ⊗ u λ + g 2 h 2 λ I d = -gh λ ∇ x B + P η g + h λ -H + λ 2 +2µ (∇ x • (hD x (u λ )) + ∇ x (h∇ x • u λ )) -κ (u λ -U ) .
(1.22) One can note the similarity between (1.22) and the classical Euler model at the low-Froude regime [START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF]. The convergence to the "incompressible" limit (λ → 0) is studied in [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF] in the context of compressible fluid flow. Besides, the approximate model (1.22) is close to the Partially Free Surface model (PFS) proposed in [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF]. Both models are able to switch between the congested areas and the free surface domain without following the interface and the surface pressure is an explicit function of the water depth. In [START_REF] Bourdarias | A finite volume scheme for a model coupling free surface and pressurised flows in pipes[END_REF], the authors link the parameter λ to a physical parameter such as the speed of sound in the water whereas it is a mathematical artifice here. As for the PFS model, one can show that the approximate model (1.22) is hyperbolic.

Proposition 1.9. Assume that ∂ t B = ∇ x B = ∂ t R = ∇ x R = κ = µ = 0, λ > 0 and h λ > 0. Then the model (1.22) is strictly hyperbolic. More precisely, denoting u N = u λ • N with a unit vector N ∈ R d-1 , the eigenvalues in the direction N are u N -c < u N < u N + c with c = 1 + 1 h λ ≥H λ 2 gh λ .
Proof. The proof is similar to the classical shallow water case. The approximate model (1.22) can be writen as ∂ t W + J x ∂ x W + J y ∂ y W = 0 with W = (h λ , u x , u y ) t , the Jacobians

J x =    u x h λ 0 1 + 1 h λ ≥H λ 2 g u x 0 0 0 u x    and J y =    u y 0 h λ 0 u y 0 1 + 1 h λ ≥H λ 2 g 0 u y    .
The eigenvalues in the direction N being the eigenvalues of the matrix J x N x +J y N y , we conclude by simple computations.

Numerical resolution

The current section is devoted to the numerical approximation of the congested shallow water model (1.8)-(1.10). Following the strategy proposed in Section 1.2, the numerical approach is based on a discretization of the hyperbolic model (1.22) at the regime λ ≪ 1. We recall that the model (1.22) came from a particular choice of L, which leads, in addition of the mandatory characteristics of the approximate model , i.e. an energy estimation, the positiveness of the relative pressure and the steady state preservation, to an explicit formula for the pressure (1.21). Another choice of L would be less efficient numerically since Newton's fixed point explained hereafter would become less straightforward. From now on, the dissipative terms are neglected, i.e. κ = µ = 0. More precisely, the numerical strategy is based on an operator splitting between the hyperbolic part of the equation and the dissipative terms. We only detail the discretization of the hyperbolic part of the model since the dissipative step might be solved using an usual implicit scheme.

The numerical strategy used to solve the approximate model (1.22) has to be consistent with the "incompressible" regime λ ≪ 1. More precisely, in the congested areas, this corresponds to high potential regime since the second term in (1.21) becomes large. By analogy with gas dynamics, this corresponds to the low-Mach regime, see [START_REF] Klainerman | Singular limits of quasilinear hyperbolic systems with large parameters and the incompressible limit of compressible fluids[END_REF][START_REF] Klainerman | Compressible and incompressible fluids[END_REF]. Classical Godunov solvers are known to be too much dissipative in this case and to require a very restrictive CFL condition, see [START_REF] Dellacherie | Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number[END_REF][START_REF] Guillard | On the behavior of upwind schemes in the low Mach number limit: II. Godunov type schemes[END_REF]. Several schemes in the literature are designed to recover the "incompressible" regime, see [START_REF] Chalons | Large time step and asymptotic preserving numerical schemes for the gas dynamics equations with source terms[END_REF][START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations[END_REF][START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF].

In addition, a particular attention is given to the preservation of the lake at rest solution. Even more important than in the case of the classical shallow water model, the conservation of the lake at rest at the numerical level, using so-called well-balanced schemes [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], is a required property for the simulation of the model (1.22).

More precisely, if a scheme is not well-balanced and does not preserve the lake at rest, it will create oscillations, of order of the space step, at the free surface. However, the scheme is still consistent with the classical shallow water model. In the case of the congested shallow water model, small oscillations can artificially separate the surface from the roof and lead to strong variations of the surface pressure.

Most of the schemes designed to be consistent in the "incompressible" regime do not care about steady states at rest. More precisely, this regime is usually analyzed for the low-Mach regime of the Euler equations without source term or with a source term constant in space such as gravity. However, the Centered-Potential Regularization (CPR) scheme introduced in [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF] satisfies the two requirements: consistency with the "incompressible" regime and the well-balanced property. The numerical strategy presented in the following is an adaptation of the CPR scheme to take into account the congestion constraint.

2.1. The CPR scheme adapted to the congested shallow water model. Let T be a mesh of Ω x composed of star-shaped control volumes, see Figure 2. We denote k ∈ T a control volume, F k the set of its faces, |k| its volume and |∂k| its surface area. Furthermore for a face f , its area is denoted |f | and k f is the neighbor control volume of k such that k ∩ k f = f . The space step is defined by l k = |k| |∂k| and the normalized face size by

ν k f = |f | |∂k| .
The unit normal to f outward to k is denoted as N k f . The set of the faces at the boundary of the mesh is designed by F ∂T . In addition, the time step is denoted by δ n+1 t , i.e. t n+1 = t n + δ n+1 t . 2.1.1. Description of the CPR scheme. The CPR scheme is an IMplicit-EXplicit Finite Volume method. The numerical unknowns are approximations of the average in a control volume k of the variables. The mean atmospheric pressure, the mean bottom level, the mean roof level and the mean opening in the control volume k and at time t n are respectively denoted by

P n k , B n k , R n k and H n k = R n k -B n k .
One particularity of the CPR scheme is that the numerical unknowns are not the conserved variables, they are chosen as the potential φ n k and the horizontal velocity u n k . This choice of variables is particularly well adapted to the congested shallow water model since the potential can be used as a parametrization of the water depth h n k and the surface pressure p n k . More precisely, we recall that the potential reads φ λ = g (h λ + B) + p λ and we set

f a c e f N k f f Cell k φ n k , u n k Cell k f φ n k f , u n k f
h (k, n, φ) =        φ -P n k g -B n k , if φ ≤ g H n k + B n k + P n k H n k + λ 2 φ-P n k g -B n k 1 + λ 2 , else (2.1)
and

p (k, n, φ) = φ -g (h (k, n, φ) + B n k ) . (2.2)
For readability, we set h n k = h (k, n, φ n k ) and p n k = p (k, n, φ n k ). Using an implicit Finite Volume method, the mass conservation, first equation of (1.22), reads

h n+1 k = h n k - δ n+1 t l k f ∈F k F n+1 f • N k f ν k f (2.3)
with F n+1 f an approximation of the mean mass flux h λ u λ at the face f between the times t n and t n+1 . More precisely if f ∈ F ∂T , we set

F n+1 f = h n+1 u n f -γδ n+1 t h n+1 l f φ n+1 k f k N k f
with a regularization parameter γ ≥ 0 characteristic of the CPR scheme, called CPR-parameter in the sequel, see [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]. We have also introduced the following notation at the face

(ψ) f = ψ k + ψ k f 2 and [ψ] k f k = ψ k f -ψ k 2 .
If f ∈ F ∂T , the numerical scheme requires boundary conditions which depend on the regime of the flow. In the current work, we do not detail the boundary conditions treatment (For a hyperbolic system it is known to be a difficult point, see [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF] for more details). The scheme (2.3) is an implicit and non-linear scheme on the potential, thus it cannot be solved directly. An iterative Newton process is used to solve it at each time step. More precisely, writing the scheme as

S k φ n+1 = 0 with S k (φ) = h (k, n + 1, φ) -h n k + δ n+1 t l k f ∈F k (h (k, n + 1, φ) u n ) f • N k f -γδ n+1 t h (k, n + 1, φ) l f [φ] k f k ν k f ,
the Newton method consists in constructing a sequence (φ n,q ) q≥0 defined by φ n,0 = φ n and J (φ n,q ) δ n,q φ = S (φ n,q ) , φ n,q+1 = φ n,q -δ n,q φ where J denotes the Jacobian matrix of S (φ) = (S k (φ)) k∈T . Furthermore φ n+1 is defined as lim q→∞ φ n,q . A similar Newton method is analyzed in [START_REF] Brenner | Improving Newton's method performance by parametrization: the case of Richards equation[END_REF]. Numerically, a quadratic convergence of the Newton method is observed.

Once the potential φ n+1 k is recovered, the water depth h n+1 k is known and the momentum balance is computed using the following explicit upwind scheme in order to get

u n+1 k h n+1 k u n+1 k = h n k u n k - δ n+1 t l k h n+1 k f ∈F k φ n+1 k f k N k f ν k f - δ n+1 t l k f ∈F k u n k F n+1 f • N k f + -u n kf F n+1 f • N k f - ν k f . (2.4)
Eventually, the discrete surface pressure can be estimated a posteriori for output using (2.2).

Remark that thanks to the definition of the variables (2.1)-(2.2), in particular the choice of φ instead of h, the scheme (2.3) and (2.4) can be used in practice with λ set to 0. In other words, there is no longer penetration in the roof and the congestion constraint (1.9) is exactly satisfied. In the next section, we first propose an analysis with λ > 0 then we investigate the case λ = 0 numerically in the simulations.

2.1.2. Numerical analysis. In the following section, the numerical counterpart of the main properties of the congested shallow water model is established. In particular, the dissipation of the mechanical energy in Proposition 1.6, acts as a mathematical entropy and proves the numerical stability of the scheme. Firstly, let us recall the main numerical properties inherited from the CPR scheme, before showing the properties specific to the congested shallow water model.

Let us start by a consistency result valid on a regular cartesian grid. We denote by δ x the uniform space step, i.e. the distance between the center of two neighboring cells. In this particular case the coordinates of the center of each cell are denoted by X k . Proposition 2.1 (Consistency). Let λ and γ be fixed. Then for smooth enough solutions and a regular cartesian grid we have

|h n k -h λ (X k , t n )| = O (δ t , δ x ) and h n k u n k -h λ (X k , t n ) u λ (X k , t n ) = O (δ t , δ x ) .
Proof. The proof is done in [40, Proposition 2.2].

Proposition 2.2. The scheme (2.3) and (2.4) satisfies the following properties i) (Positivity of water depth) Assume h 0 k > 0 for any k ∈ T and assume that the following implicit CFL condition is satisfied for any n ∈ N

u n f • N k f + γ 2 [φ n+1 ] k f k δ n+1 t ≤ min h n+1 k , h n+1 k f h n+1 k + h n+1 k f min l k , l k f . (2.

5)

Then h n k > 0 for any k ∈ T and n ∈ N. ii) (Conservation of water volume) For any n ∈ N,

k∈T |k| h n+1 k = k∈T |k| h n k -δ n+1 t f ∈F ∂T F n+1 f • N k f |f |.
iii) (Well-balanced scheme) For any φ 0 ∈ R, let the initial condition be φ 0 k = φ 0 and u 0 k = 0. Then for any time, the solution of (2.3) and (2. One can check that the discrete lake at rest iii) defined in Proposition 2.2 is a discrete version of the lake at rest of Proposition 1.8.

Note that the CFL condition (2.5) cannot be satisfied at a wet-dry interface. In practice the simulation of the latter can be realised by taking artificially some small time step but the scheme will not satisfy the dissipation law given below in Proposition 2.3.

Let us now show some properties in the case of the congested shallow water model.

Proposition 2.3. For any γ ≥ 1 and under the CFL condition (2.5), the CPR scheme (2.3) and (2.4) satisfies the following dissipation law

E n+1 k + δ n+1 t l k f ∈F k G n+1 K,f + G n+1 P,f ν k f ≤ E n k + S n+1 k
with the discrete mechanical energy E n k = K n k + P n k , the kinetic energy and the potential energy respectively given by

K n k = 1 2 h n k u n k 2 and P n k = gh n k B n k + h n k 2 + h n k P n k + g 2λ 2 h n k -H n k 2 +
, the flux of the kinetic energy

G n+1 K,f = 1 2 u n k 2 F n+1 f • N k f + - 1 2 u n k f 2 F n+1 f • N k f - + δ n+1 t 2 h n+1 k f l k f - h n+1 k l k φ n+1 k f k 2
, the flux of the potential energy

G n+1 P,f = φ n+1 f F n+1 f • N k f + h n+1 k u n k -h n+1 k f u n k f 2 • φ n+1 k f k N k f k
and the source term

S n+1 k = -p n+1 k -P n+1 k H n+1 k -H n k + h n k gB n+1 k + P n+1 k -gB n k -P n k .
Proof. The following proof is an adaptation of [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]Theorem 2.3]. The estimation of the discrete kinetic energy dissipation law under the CFL condition (2.5) established in [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]Lemma 2.4] still holds in our case since it is obtained for any potential. However, the discrete potential energy dissipation law [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]Lemma 2.5] is obtained for a differentiable potential, which is not our case due to the positive part function of the regularization. Let us establish the discrete potential energy balance in our case. Using the two equalities b

+ (b -a) = b + (b + -a + ) + b + a -where a -= a + -a and 2a (a -b) = a 2 -b 2 + |a -b| 2 , simple calculations lead to φ n+1 k h n+1 k -h n k = P n+1 k -P n k + Q n+1 k -S n+1 k where Q n+1 k = g 2 h n+1 k -h n k 2 + g 2λ 2 h n+1 k -H n+1 k + -h n k -H n k + 2 + g λ 2 h n+1 k -H n+1 k + h n k -H n k - . Now multiplying (2.3) by φ n+1 k
and taking into account that Q n+1 k ≥ 0 gives the discrete potential energy dissipation law

P n+1 k ≤ P n k - δ n+1 t l k f ∈F k G n+1 P,f -R n+1 P,f ν k f + δ n+1 t W n+1 k + S n+1 k
with the rest coming from the space discretization and the work of the conservative forces respectively given by

R n+1 P,f = -γδ n+1 t h n+1 l f φ n+1 k f k 2 and W n+1 k = h n+1 k u n k l k • f ∈F k φ n+1 k f k N k f ν k f .
Finally, the proof of [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]Theorem 2.3] allows to conclude.

Note that S n+1 k δ n+1 t is consistent with the terms depending on a time derivative in the right-hand side of the energy balance law given in Proposition 1.6.

Proposition 2.4. The scheme (2.3) and (2.4) satisfies the following properties i) (Positivity of relative pressure) The relative pressure is positive, i.e. for any k ∈ T and n ∈ N

1 -1 h n k <H n k (p n k -P n k ) = 0 and p n k ≥ P n k .
ii) (Upper bound on the surface pressure) Assume that for any f ∈ F ∂T , we have F n f = 0. Then there exists K ∈ R * + such that for any λ > 0, we have p n k ≤ K for any k ∈ T and n ∈ N.

Proof. The proof of i) is a consequence of Proposition 1.7.

Let us now prove ii). The mass conservation (2.3) at time t n+1 reads

-γδ n+1 t f ∈F k h n+1 l f φ n+1 k f k ν k f + l k h n+1 k δ n+1 t = l k h n k δ n+1 t - f ∈F k h n+1 u n f • N k f ν k f .
Let us consider the latter equality for all k ∈ T and denote φ n+1 the vector with the components φ n+1 k for k ∈ T. We get a system of the form M n+1 φ n+1 = Φ with M n+1 a non-singular M-matrix. Indeed the entries of the main diagonal of M n+1 are non-negative, the off-diagonal entries are non-positive and the matrix is strictly diagonally dominant. The system leads to an inequality of the form

M n+1 φ n+1 ≤ Ξ n+1 with Ξ n+1 = l k h ∞ δ n+1 t + h ∞ u n k .
In fact, thanks to the water volume conservation ii) of Proposition 2.2, the water depth is bounded by 0

≤ h n k ≤ h ∞ with h ∞ = k∈T |k|h 0 k min k∈T |k|
and the bound on the right hand side follows.

As a consequence, we get a bound on φ n k and the bound on the surface pressure follows. More precisely, ∀n ≥ 0, max

k∈T p n k ≤ max k∈T φ n k + g (h ∞ + max k∈T B n k ).
Thanks to ii) from Proposition 2.4, we are able to characterize the penetration with respect to λ. More precisely considering (1.21), we have

h n k -H n k + = O λ 2 .
It shows that (1.9) is verified as λ goes to zero. Furthermore, it shows that the mechanical energy in Proposition 2.3 remains bounded when λ tends to 0. Unless indicated otherwise, for all the test cases g = 9.81, P η = 0, γ = 1 and a channel of length l = 1 is considered, i.e. x ∈ [0, 1] where the entry is located at x = 0 and the exit at x = 1.

2.2.1.

The lake at rest. The discrete lake at rest iii) of Proposition 2.2 is tested and validated in various configurations with bottom and roof which may be continuous or not. Note that in such configurations, the CFL condition (2.5) is satisfied for any positive time step. Setting λ to zero, the discrete lake at rest iii) of Proposition 2.2 satisfies the constraint (1.9) and thus is a discrete version of Proposition (1.14). This setting is tested and validated as well. 2.2.2. Transcritical steady flow with free hydraulic jump. For any smooth enough bottom and roof, an analytical solution for transcritical steady flow can be obtained as long as the top of the hydraulic jump does not reach the roof. The description of the analytical solution and the algorithm used to compute the solution are described in Appendix A.

Let us consider a roof and a bottom defined by R (x) = 100 (x -0.5) 2 + 5.8 and

B (x) = 1, if x ≤ 0.65 exp -100 (x -0.65) 2 , if x > 0.65.
The water depth at the exit is fixed to h (l, t) = 8 and the discharge at the entry is set to hu (0, t) = 40. The physical parameters are chosen such that the flow is subcritical (fluvial) at the boundaries but becomes supercritical (torrential) because of the roof. The numerical solutions with the parameter λ = 0 and λ 2 = δ x together with the analytical solution with regular space step δ x = 10 -3 are plotted in Figure 3.

First of all, let us discuss the results with respect to the value of λ. The simulation was performed with several parameters λ from 1 to 0 and for δ x ∈ { 1 /300, 1 /1000, 1 /3000, 1 /10000}. In Figure 4 the L 2 -errors of the water height, the velocity and the surface pressure are plotted versus λ 2 for the different space steps. We notice that for λ < 1, the water depth and the velocity are well approximated. The surface pressure is more difficult to obtain. The best agreement with the analytical solution is obtained for

δ x ≤ λ 2 ≤ (δ x ) 1 2
. For λ 2 < δ x , the surface pressure is not in good agreement with the analytical solution. Indeed the numerical error of the scheme on the congestion constraint is fully reflected to the surface pressure when λ tends to 0 whereas it is split between the water level and the surface pressure in the case λ > 0.

Figure 5 shows the convergence rates in L 2 -norm to the analytical solution for the water depth, the velocity as well as the surface pressure for λ = 0 and λ 2 = δ x . In both cases, λ = 0 and λ 2 = δ x , the water depth and the velocity converge at first order. In the case λ 2 = δ x , the surface pressure converges at approximatively order 1.5, which is faster than expected. Indeed, as the water height converges at first order, the surface pressure is expected to converge at first order as well. In the case λ = 0, the surface pressure seems to converge but chaotically.

2.2.3.

Transcritical steady flow with constrained hydraulic jump. The solution presented in Appendix A is not valid when the top of the hydraulic jump reaches the roof. In this case, the non-conservative product h∇ x p η is not clearly defined. The following simulation illustrates the solution of the numerical scheme in such a configuration.

Let us consider a bottom and a roof respectively defined by B (x) = 3.5 exp -100 (x -0.35) 2 and R (x) = 100 (x -0.5) 2 + 5.8.

The water depth at the exit is fixed to h (l, t) = 8 and the discharge at the entry is set to hu (0, t) = 25. The physical parameters are chosen such that the flow is subcritical (fluvial) at the boundaries but becomes supercritical (torrential) because of the nonflat bottom and the hydraulic jump reaches the roof.

The simulated solution with δ x = 10 -3 and λ 2 = δ x as well as λ = 0 is plotted in Figure 6. Furthermore a reference solution is plotted for comparison. More precisely the analytical solution of the water height without considering the roof and the corresponding velocity are plotted, see (A.2) as well as the analytical value of the surface pressure considering the roof, see (A.3). As expected, the water depth and the velocity are similar in the areas where the flow does not reach the roof and only the localization of the hydraulic jump is not recovered. At the hydraulic jump, the water depth as well as the surface pressure are discontinuous which makes the product h∇ x p η really non-conservative and the position of the shock not easy to compute. The peak of the pressure at the top of the hydraulic jump is due to oscillation effects of the CPR scheme at discontinuities, see [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]. By increasing the CPR-parameter γ in order to regularize the solution and decrease the oscillation effects, the approximation of the water height gets worse. However, the hydraulic jump is clearly not located at the same position with and without the roof. We conclude that the surface pressure has an impact on the jump condition, at least at the discrete level. 2.2.4. Time-dependent emptying tank. Assuming a fixed and flat bottom, an analytical solution of the congested shallow water model can be obtained for any fixed bottom as long as the interface between the congested and the free surface flow is fixed. More precisely, we claim the following statement.

Proposition 2.5. Let us consider a fixed flat bottom and an uniform atmospheric pressure, i.e. B = ∂ x P = 0. Furthermore, let us consider a fixed constant roof in a part of the domain delimited by two fixed points x l ≤ 0 ≤ x r , i.e.

R (x) = H 0 , if x l ≤ x ≤ x r ∞,
elsewhere with an arbitrary opening H 0 > 0.

Then for any α ≥ 0, h 0 > H 0 and any time t ≤ h 0 -H0 αH0 , the following functions Proof. The proof relies on simple computations after injecting the functions in the congested shallow water equations (1.8)-(1.10). The non-conservative product is well-defined since the pressure at the bottom gh + p η is continuous.

h (x, t) =    H 0 , if x l ≤ x ≤ x r h 0 1 + αt , elsewhere, u (x, t) =          α (x -x g ) 1 + αt , if x < x l 0, if x l ≤ x ≤ x r α (x -x d ) 1 + αt , if x r < x p η (x, t) = P η (t) + g h 0 1 + αt -H 0 0.
Figure 7 shows the analytical solution of the Proposition 2.5 at time t = 0.5 with the parameters h 0 = 10, H 0 = 5.8, x l = 0.3 and x r = 0.7, together with the numerical solution with δ x = 10 -3 and λ 2 = δ x as well as λ = 0. The numerical solution is in very good agreement with the analytical solution. Figure 8 shows the convergence rates in L 2 -norm to the analytical solution for the water depth, the velocity and the surface pressure for λ = 0 and λ 2 = δ x . The three unknowns converge approximately at order 1.36 for any λ, which is better than the expected first order. In this case the pressure depends linearly on the water depth, which explains the same order of convergence for the pressure for both values of λ compared to the stationary case. 2.2.5. Time-dependent filling tank. The current simulation is a numerical investigation of the jump condition when the top of the hydraulic jump reaches the roof. More precisely, a fixed flat bottom and roof and a symmetrical initial condition given by

h 0 (x) = 1 2 H (x, 0) and u 0 (x) = C, if x < 1 2 -C, if x > 1 2
with H (x, 0) = 5 for all x ∈ [0, l] and C = 8 are considered. The initial condition is chosen such that the flow is supercritical (torrential). Figure 9 shows the numerical result at time t = 0.5 with δ x = 10 -4 , λ 2 = δ x and γ = 10. Due to the shock waves, the CPR-parameter γ is increased to get slightly more regularisation, see [START_REF] Parisot | Centered-potential regularization for the advection upstream splitting method[END_REF]. Otherwise oscillations of the pressure in the vicinity of the shock due to oscillation effects of the CPR scheme are observed. However, the oscillations increase when λ becomes smaller and more regularization is needed. As a consequence, we cannot pass λ = 0 in the simulation, as either the oscillations become too large or the regularization effect is too high. As expected, the solution is composed by tree constant states separated by two shock waves. By symmetry, it is clear that the velocity in the middle state should vanish. In addition, the conservation of the mass imposes the velocity of the shock wave. More precisely, it leads to X ′ = 2U with X (t) the position of the shock. However, the computation of the pressure of the middle state is inaccessible since the water depth and the surface pressure are discontinuous at the shock.

Conclusion

A shallow water type model with a constraint of congestion is derived from the Navier-Stokes equations for the modeling of a roof in geophysical applications. The simulation of the congested shallow water model is challenging due to the unknown localization of the interface between the congested areas and the free surface domain. To overcome this difficulty, an approximate model which is hyperbolic in the free surface domain as well as in the congested areas is proposed based on a pseudocompressible approximation. We show that two properties, mainly well-balanced and low-Froude stability, are needed at the discrete level for the simulation of the approximate model. Eventually, a numerical scheme adapted from the Centered-Potential Regularization scheme is proposed and confronted with analytical solutions to illustrate the efficiency of the scheme.

Concerning the theoretical analysis, several issues have not yet been adressed. Firstly, to confirm the relevance of the pseudo-compressible approximation, the convergence of the approximate model (1.16) and (1.18) as λ goes to zero to the initial congested shallow water model (1.8)-(1.10) seems a natural question but not at hand for the moment. The second open problem in the analysis of the congested shallow water model (1.8)-(1.10) is the definition of the non-conservative product h∇ x p η . More precisely the water height as well as the surface pressure can be discontinuous in the same point as shown in Figure 6. Since the surface pressure plays the role of a Lagrange multiplier, the model (1.8)-(1.10) is not a hyperbolic conservation law. As a consequence a classical analysis of non-conservative products in hyperbolic theory, for example [START_REF] Dal Maso | Definition and weak stability of nonconservative products[END_REF], cannot be directly applied here.

To improve the modeling, several ways can be investigated. First of all, the modeling of areas with negative relative pressure, i.e. p η < P η is quite challenging. This issue seems more generally linked to the modeling of trapped air pockets under the roof. Secondly, the congestion constraint can be applied to more complex models of free surface flow such as the layerwise discretized models [START_REF] Audusse | A multilayer Saint-Venant model: derivation and numerical validation[END_REF][START_REF] Audusse | A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation[END_REF], the non-hydrostatic models [START_REF] Bristeau | An energy-consistent depth-averaged euler system: derivation and properties[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF] or the combination of the two [START_REF] Fernandez-Nieto | A hierarchy of non-hydrostatic layer-averaged approximation of Euler equations for free surface flows[END_REF]. Moreover, the bilateral interaction between the roof and the fluid seems interesting for applications, in particular for object buoyancy such as boats or buoys. Last but not least, the modeling of an immersed object, using layerwise discretized models, is still quite far but can be a challenging objective.

Appendix A. The transcritical analytical solution

The classical test case of transcritical steady flow is largely used to measure the convergence of a numerical scheme confronted to a hydraulic jump. Even if an explicit formula is not obtained, the exact solution can be computed on a grid. The analytical solution for free surface flow is described in [START_REF] Chanson | Hydraulics of Open Channel Flow[END_REF]. This section is devoted to the analytical solution of the congested shallow water model. However, since the non-conservative product is not clearly defined, the strategy is only valid assuming the solution is continuous when the roof is reached. In other words, the top of the hydraulic jump is supposed not to reach the roof. In addition, the boundary conditions are given for the classical shallow water model in the subcritical regime. It follows that the flow does not reach the roof at the boundary.

The solution is given by intervals with constant hydraulic heads separated by discontinuities satisfying the classical Rankine-Hugoniot relation since the hydraulic jump is supposed not to reach the roof. More precisely, for steady flow, the mass conservation equation of (1.8) leads to a constant mass flux hu = Q. Moreover, when the water depth is continuous, the energy conservation Proposition 1.3 holds and yields a constant energy flux 1 2 u 2 + g (h + B) + p η hu = G. Assuming that the surface pressure and the water depth are known in a point, denoted from now on characteristic point X, the hydraulic head K = G Q can be computed. More precisely, the hydraulic head reads

K = h (X) + B (X) + p (X) g + H 3 c 2h 2 (X) (A.1)
with

H c = Q √ g 2 3
the critical water depth.

At the free surface, the surface pressure corresponds to the atmospheric pressure and the water depth is solution of the following third order polynomial p η = P η and

h 3 + B + P η g -K h 2 + H 3 c 2 = 0. (A.2)
More precisely, one can show that the water depth is the larger root of the polynomial in the area where the flow is subcritical and the smaller of the positive roots when the flow is supercritical. In the same interval, if the roof is reached, the water depth corresponds to the opening and the surface pressure is defined by the following expression The surface pressure in (A.3) is a generalization of the Archimedes' principle since when the mass flux vanishes, i.e. Q = 0, thus H c = 0, the surface pressure corresponds to the weight of a fictitious water volume above the roof. Now let us explain the strategy followed to compute the analytical solution on a grid. In the following, we only detail how to determine the characteristic node j, discrete counterpart of the characteristic point X.

p η = g K -H -B - H 3
In the free surface area, the characteristic node j is determined by the point where the water depth corresponds to the critical water depth. Indeed, in this point the water depth is known (the pressure is determined by the atmospheric pressure) and the hydraulic head K can be computed using (A.1). More precisely, as it is explained in [START_REF] Chanson | Hydraulics of Open Channel Flow[END_REF], this point corresponds to the point where B + Pη g is locally decreasing.

In the congested area, the characteristic node j is determined by the right point where the flow reaches the roof. Indeed since this point reaches the roof, we have h j = H j and since at its right the flow has a free surface and the surface pressure is assumed continuous we have p j = P j . To determine this point, we start from the right point where the opening intersects the critical water depth. At the left of this point, the flow is clearly supercritical. Then, we proceed successively node by node going to the left as follows. Assuming the current node is the characteristic node, the hydraulic head is computed using (A.1). Then, the pressure at the left neighbor node is computed using (A.3). If the relative pressure at the left is negative, we conclude that the current node is not the characteristic node and we pass to the next one. Once the characteristic node is determined, the water depth at the right is computed using (A.2) until the hydraulic jump. At the left of the characteristic node, the surface pressure is computed using (A.3) until the pressure becomes negative, which corresponds to the entrance of the water under the roof. Eventually, the velocity can be deduced dividing in each point the mass flux by the water depth. The numerical implementation of the analytical solution is given in Algorithm 1.

1 Figure 1 .

 11 Figure 1. Geometrical description of the flow

1. 1 . 1 .

 11 Navier-Stokes model with roof. A Navier-Stokes model describing the dynamics of the fluid is proposed.

Proposition 1 . 3 .

 13 Any smooth enough solution of the congested shallow water model (1.8)-(1.10) satisfies the following energy balance law

Figure 2 .

 2 Figure 2. Numerical notations for the two dimensional case, i.e. d = 3

  4) reads φ n k = φ 0 and u n k = 0. Proof. The proofs of i), ii) and iii) are done in [40, Proposition 2.1, Proposition 2.6, Proposition 2.8].

Figure 3 .

 3 Figure 3. Transcritical steady flow with free hydraulic jump.

Figure 4 .L 2 -Figure 5 .Figure 6 .

 4256 Figure 4. Transcritical steady flow with free hydraulic jump: L 2errors for the water height, the velocity and the surface pressure versus the parameter λ for different values of the space step δ x

Figure 7 .L 2 -Figure 8 .Figure 9 .

 7289 Figure 7. Emptying tank at time t=0.5

  Sections 2.2.1, 2.2.2 and 2.2.3 are devoted to stationary solutions. More precisely, simulations of the lake at rest are commented in Section 2.2.1. In Section 2.2.2, simulations of transcritical steady flow with a free hydraulic jump are compared to an analytical solution. A simulation of a transcritical steady flow with an hydraulic jump constrained by the roof is illustrated in Section 2.2.3. Then Sections 2.2.4 and 2.2.5 are devoted to time-dependent solutions. In Section 2.2.4, a time-dependent analytical solution of (1.8)-(1.10) is given and compared to the numerical solution obtained with (2.3)-(2.4). Eventually, a simulation in the case of discontinuous velocity is realized in Section 2.2.5.

2.2. Simulation and numerical validation. The numerical solution of the congested shallow water model is illustrated in a one dimensional framework, i.e. x ∈ R.

  N -integer scalar : Number of control volumes H c -real scalar : Critical water depth H r -real scalar : Right water depth (boundary condition) B -real array(1:N) : Bottom H -real array(1:N) : Opening P -real array(1:N) : Atmospheric pressure g -real scalar : Gravity accelerationfor i=1 to N do (h i , p i ) ← (H c , P i ); # Initialization # Starting from the right border (h N , p N ) ← (H r , P N ) ; K ← (A.1) ; subcritical ← (h N > H c ) ; i ← N -1 ; while (i > 0) do if subcritical then (h i , p i ) ← (A.2) ; # Assuming free surface if h i > H i then (h i , p i ) ← (A.3) ; # If congested else # Find the characteristic point j ← i + 1 ; f ound ← f alse ; while not f ound do j ← j -1 ; if (j = 0) then STOP 42 # Left boundary supercritical if H j ≥ H c then if (gB j-1 + P j-1 ≤ gB j + P j ) then (h j , p j ) ← (H c , P RH ← f alse ; while not RH do i ← i + 1 ; (h i , p i ) ← (A.2) ; # Assuming free surface if h i > H i then (h i , p i ) ← (A.3) ; # If congested if (i = N ) then RH ← true ; # Right boundary supercritical -1 ; (h i , p i ) ← (A.3); while (p i > 0) do i ← i -1 ; (h i , p i ) ← (A.3) ; subcritical ← (h i > H c ) ; (h i , p i ) ← (A.2) ; if h i < H c then (h i , P i ) ← (H c , P i ); subcritical ← (h i > H c ) ; i ← i -1 ;# Test the regime for i=1 to N do u i ← √

	Algorithm 1: Transcritical analytical solution
	Input:
	else RH ←	H 3 c hi +	h 2 i 2 <	H 3 c hi+1 +	h 2 i+1 2	;	# Rankine-Hugoniot
	#	Left congested side until free surface
	i ← j gH 3 c hi ;				# Computation of the velocity
				c 2H 2		and	h = H.	(A.3)

j ) ; K ← (A.1) ; f ound ← true ;

# Found free surface area else (h j , p j ) ← H j , P j ; K ← (A.1) ; (h j-1 , p j-1 ) ← (A.3) ; if (p j-1 > P j-1 ) then f ound ← true ; # Found congested area # Right side until free hydraulic jump i ← j -1 ;

Acknowledgments

The authors wish to express their warm thanks to Cindy Guichard for many fruitful discussions and kind advices. This research is part of a collaborative program with the ANR-FEM (France Energies Marines) project HYFLOEFLU through the project call ITE EMR 2015.


Daniel Mariuta

Lucien Baldas
email: lucien.baldas@insa-toulouse.fr

Jürgen Brander

Stéphane Le Calvé

Stéphane Colin

Pascale Magaud

Christine Barrot-Lattes

A Marie-Curie-ITN within H2020 MIGRATE2017: 154722 GAS-LIQUID MICRO-SEPARATORS FOR TRAPPING OF VOLATILE ORGANIC COMPOUNDS

Keywords: Microfluidics, two-phase flow, mass transfer, uptake yield, real time detection, ultra-portability, miniaturization, high sensitivity, formaldehyde, air-borne pollution, analyzer

ABSTRACT

Volatile Organic Compounds (VOCs) are substances which were proved to be harmful for the human health. People are spending about 90% of their time indoors where these chemicals are present in high concentrations [1]. Among VOCs, formaldehyde is a major pollutant due to multiple indoor pollution sources (furniture, paintings, etc.) and high degree of toxicity being classified internationally as carcinogenic. Starting from the necessity of real time measurement of formaldehyde concentration, a device has already been developed by both ICPEES and In'AirSolutions in Strasbourg [2] which was proved to be sensitive and portable.

The analyzer works following three main stages. In the first stage, an annular two-phase flow is established in order to trap gaseous formaldehyde into an aqueous solution containing a chemical reagent. In the second stage, the kinetics of the reaction between formaldehyde and reagent are accelerated inside an oven at a specific temperature. Resulted chemical reaction product is a fluorescent compound and its concentration is measured inside a fluorescent cell, in the third stage of the process. Separation techniques using two-phase flows in micro-and mini-channels turned to be much more interesting since the increasing scientific interest in the fields of microfluidics and microfabrication techniques registered a consistent advancement. Among all applications of microfluidic separation, utilization of this technique for identification of air pollutant concentrations is very important due to its high miniaturization potential. For being capable of promoting merchandisable micro-separation systems, it is necessary to develop microsystems working in stabilized and reproducible conditions. This work is part of a common effort for improving and miniaturizing the capabilities of this existing formaldehyde analyzer to integrate it into a real time miniaturized, ultra-portable, autonomous and reliable formaldehyde measurement system. More precisely, two areas for improvement have been identified. On one hand, a new approach for improving the uptake yield of the gaseous formaldehyde has to be proposed in order to diminish the response time and the liquid reagent consumption. On the other hand, the current fluorescent detection cell comes with a large dead volume and it is necessary to develop another design configuration as part of the improvement of the detection process. Some commercial fluorescent cells with different volumes have been previously tested and different problems affecting either the response time or the quality of the measured formaldehyde concentration have been identified.

In order to identify new solutions for improving the uptake yield, a detailed knowledge of the flow conditions is required. Despite all the progress registered until now, not all the phenomena regarding two phase micro-flows are completely understood, many aspects being under investigation and closely interconnected to the progress registered in adjacent fields, such as flow visualization and numerical techniques [START_REF] Sobieszuk | Investigations of mass transfer in annular gas-liquid flow in microreactor[END_REF]. Flow patterns are strongly affected, indeed, by a multitude of parameters at microscale, e.g. surface roughness, material properties, geometry of the inlet, channel diameter, channel cross section, superficial velocities, fluids properties [START_REF] Chinnov | Regimes of two-phase flows in micro-and minichannels[END_REF]. Among all two-phase flow patterns, slug and annular patterns are utilized in the majority of the proposed separation processes in the literature due to large interfacial area and high uptake yield. In order to establish which is the suitable compromise between all the characteristics proposed by the two patterns, they have to be experimentally analyzed. The interplay of all effects involved have to be controlled and fully understood in order to propose a highly efficient and reproducible micro-separation method.

Regarding the detection stage, we are interested in developing and testing a new fluorescent integrable detection configuration in order to avoid large liquid dead volumes. It should respect some criteria such as low-cost, low detection limit, long term stability and working properly in a broad range of environmental conditions. There are some attempts in the literature to develop integrated fluorescence micro-detectors for in time measurement [START_REF] Chabinyc | An integrated fluorescence detection system in poly (dimethylsiloxane) for microfluidic applications[END_REF][START_REF] Wang | A portable time-domain LED fluorimeter for nanosecond fluorescence lifetime measurements[END_REF], fluorescence being largely used in lifetime measurements due to its selectivity, sensitivity and multiplexing capabilities. A prototype is under development in this moment in our team and the challenge consists in overcoming the constraints related to the microfabrication limitations, inlet configuration (which many times is a source of bubble entry), light pathways, LED/photodiode miniaturization and materials used for the cell. Polymeric materials and especially polydimethylsiloxane comes with promising characteristics regarding the possibility of developing embedded fluorescent systems for continuous monitoring. PDMS is optically transparent over the entire visible range and over a part of the UV range, permits rapid prototyping and scalable manufacturing with lower costs comparing the other materials used in microfluidics such as glass or silicone, which are expensive.

A preliminary experimental set-up is under development in order to visualize and characterize the specific parameters of the two-phase flow. Accordingly with the state of the art in the field of two-phase microflows, we are interested in quantifying parameters which are relevant and influential for the uptake process such as gas and liquid flow rates, gas-liquid-solid interface behavior, temperature, material properties (wettability), pressure drop. Temperature influence is very important due to our interest in coupling the first two stages of the analyzer. Heat transfer effect will be studied and its influence over the uptake yield will be emphasized and carefully analyzed. These steps will provide us with important information being helpful in the process of evaluation of mixing, mass and heat transfer and reaction yield. Iteration and evaluation of the reaction zone by quantitative criterion and correlations is a key method in designing methodology of microintegrated structures.

Experimental campaigns will be conducted for different experimental conditions (e.g. different materials, diameters and cross sections for the micro-channels, different flow patterns and flow rates, heating coupled with the uptake stage) in order to study the uptake yield. As a first step, the experiments will be conducted using fused silica micro-channels with different diameters ranging from 50µm to 530µm in order to quantify the results through empirical and optical methods. Starting from the experimental set-up and results proposed by [START_REF] Guglielmino | On-line gaseous formaldehyde detection by a microfluidic analytical method based on simultaneous uptake and derivatization in a temperature controlled annular flow[END_REF] regarding the efficiency of the uptake yield which is increased with the interior diameter, we are planning to characterize the flow dynamics in these circumstances in order to be able to propose an improved trapping system. The possibility of integrating experimental colorimetric methods, as for example the one proposed by Yang et. al. [START_REF] Yang | Visualization and characterization of gas liquid mass transfer around a Taylor bubble right after the formation stage in microreactors[END_REF], will be analyzed for a better visualization of mixing hydrodynamics and interface interaction. Further, according to the conclusions of the experimental campaigns, numerical analysis using Direct Simulation Monte Carlo (DSMC) will be conducted to characterize diffusion kinetics at molecular level. In parallel, according with the results of the bibliographical review about fluorescent micro-integrated systems, a new prototype will be proposed, analyzed and iteratively improved in order to acquire the desired capabilities.

Acknowledgements

This project has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Sklodowska-Curie Grant Agreement No. 643095.