Raphaël Cerf

Carlo Mariconda
email: mariconda@pdmat1.math.unipd.it

ORIENTED MEASURES

Keywords:

įs positive. We study the range R of an oriented measure:

: χ E has less than n -1 discontinuity points } .

1991 Mathematics Subject Classification

However for a given ρ, 0 ≤ ρ ≤ 1, the usual proofs based on convexity-extreme points arguments [START_REF] Halmos | The range of a vector measure[END_REF][START_REF] Olech | The Lyapunov theorem: its extensions and applications[END_REF] do not give any information about the existence of a "nice" set E such that

µ(E) = b a ρ dµ.
Consider for instance the two-dimensional vector measure µ(A) = (|A|, |A| + 2|A ∩ B|) where B is a borelian subset of [a, b] and | | denotes the Lebesgue measure. For each set E, the equality µ(E) = µ(B) implies B = E. When the measure µ admits a density f , Halkin [START_REF] Halkin | On a generalization of a theorem of Lyapounov[END_REF] showed that if for each vector p ∈ R n the set t ∈ [a, b] : p • f (t) > 0 (where • is the usual scalar product) is a finite (respectively countable) union of intervals then there exists a set E which is a finite (resp. countable) union of intervals. In our paper [START_REF] Cerf | Oriented measures with continuous densities and the bang-bang principle[END_REF] we introduced the stronger orientation condition ∆: we say that n real functions f 1 , • • • , f n verify condition ∆ on an interval [a, b] if for each k in {1, • • • , n}, the determinant

f 1 (x 1 ) f 1 (x 2 ) • • • f 1 (x k ) f 2 (x 1 ) f 2 (x 2 ) • • • f 2 (x k ) . . . . . . . . . . . . f k (x 1 ) f k (x 2 ) • • • f k (x k )
is not equal to zero whenever the x i 's in [a, b] are distinct and its sign is constant on the k-tuples (x 1 ,

• • • , x k ) such that a ≤ x 1 < x 2 < • • • < x k ≤ b.
We showed that if a measure µ admits a density function whose components are continuous and satisfy the orientation condition ∆ then the set E may be built in such a way that its characteristic function has at most n discontinuity points. Moreover, if 0 < ρ < 1 there exist two such sets E 1 and E 2 whose characteristic functions χ E 1 and χ E 2 have exactly n discontinuity points (one set is a neighbourhood of a whereas the other is not). Our proofs relied upon the fact that the map

(α 1 , • • • , α n ) -→ α 2 α 1 f (x) dx + α 4 α 3 f (x) dx + • • •
is differentiable and has an invertible Jacobian whenever a < α 1 < • • • < α n < b.

We also showed that whenever a function x satisfies x(0) = • • • = x (n-2) (0) = 0 and x (n-1) (0) = 1 then the n functions (x (n-1) , • • • , x , x) verify ∆ on a neighbourhood of 0.

We applied these results to the study of reachable sets of constrained bang-bang solutions and to non-convex problems of the calculus of variations.

In this paper we deal with measures which are not necessarily absolutely continuous with respect to the Lebesgue measure.

Oriented measure. If A 1 , • • • , A k are k measurable sets of [a, b], by A 1 < • • • < A k we mean that for all k-tuple (x 1 , • • • , x k ) of A 1 × • • • × A k we have x 1 < • • • < x k . A measure µ = (µ 1 , • • • , µ n ) is said to be oriented if for each k-tuple of measurable sets A 1 , • • • , A k such that A 1 < • • • < A k the determinant µ 1 (A 1 ) • • • µ 1 (A k ) . . . . . . . . . µ k (A 1 ) • • • µ k (A k ) is positive.
In this more general framework we give a new proof of the results stated in [START_REF] Cerf | Oriented measures with continuous densities and the bang-bang principle[END_REF].

We carry out a deep study of the range R of the measure:

• for each point q of its interior R • there exist exactly two distinct "dual" sets E 1 , E 2 whose characteristic functions have n discontinuity points such that µ(E 1 ) = q = µ(E 2 );

• the set R • coincides with b a ρ dµ : 0 < ρ < 1 so that the above set is open;

• the set R is strictly convex;

• a point µ(E) belongs to the boundary ∂R of R if and only if the characteristic function of E has less than n -1 discontinuity points;

• finally we give a recursive decomposition of the boundary ∂R.

Oriented measures

Throughout the paper we will work with an interval [a, b] equipped with the Lebesgue σfield L. Measurable will mean measurable with respect to this σ-field. A negligible set is a measurable set of Lebesgue measure zero. A vector measure on [a, b] is a countably additive set function defined on the Lebesgue σ-field with values in R n for some integer n.

Notation. If A 1 , • • • , A k are k measurable sets of [a, b], by A 1 < • • • < A k we mean that A 1 , • • • , A k have non-zero Lebesgue measure and for all k-tuple (x 1 , • • • , x k ) of A 1 ו • •×A k we have x 1 < • • • < x k . Let µ = (µ 1 , • • • , µ k ) be a vector measure. If ρ belongs to L 1 µ ([a, b]), we note µ i (ρ) = b a ρ dµ i , µ(ρ) = b a ρ dµ = µ 1 (ρ), • • • , µ k (ρ) . Definition 2.1. A vector measure µ = (µ 1 , • • • , µ n ) on [a, b] is said to be oriented on [a, b] if it is non-atomic and if for each k in {1, • • • , n} and for each k-tuple of measurable sets A 1 , • • • , A k such that A 1 < • • • < A k the determinant µ 1 (A 1 ) • • • µ 1 (A k ) . . . . . . . . . µ k (A 1 ) • • • µ k (A k ) is positive.
Remark. If µ is oriented then µ 1 is a positive measure which assigns positive values to sets of positive Lebesgue measure. In particular, if I is a non-trivial interval, then µ(I) is non-zero.

Remark. If µ is oriented and

I 1 , • • • , I n are n disjoint non-trivial intervals, then the vectors µ(I 1 ), • • • , µ(I n ) form a basis of R n .
A very important fact concerning oriented measures is that their characteristic property carries on from sets to positive functions.

Notation. If ρ is a function its support is the set supp ρ = { x : ρ(x) = 0 }.

Theorem 2.2. Let µ = (µ 1 , • • • , µ n ) be an oriented measure. If ρ 1 , • • • , ρ n are n µ integrable non-negative functions such that supp ρ 1 < • • • < supp ρ n then the determinant µ 1 (ρ 1 ) • • • µ 1 (ρ n ) . . . . . . . . . µ n (ρ 1 ) • • • µ n (ρ n ) is positive.
Let us first state a preparatory lemma.

Lemma 2.3. Let µ = (µ 1 , • • • , µ n ) be a vector measure and ρ 1 , • • • , ρ n be n measurable µintegrable functions. The determinant

ρ 1 dµ 1 • • • ρ n dµ 1 . . . . . . . . . ρ 1 dµ n • • • ρ n dµ n is equal to • • • ρ 1 (s 1 ) • • • ρ n (s n ) d σ∈S n (σ)µ σ(1) ⊗ • • • ⊗ µ σ(n) (s 1 , • • • , s n ).
Proof of the lemma. The identity is obviously true whenever ρ 1 , • • • , ρ n are characteristic functions. The monotone class theorem yields the result.

Proof of theorem 2.2. We apply the lemma. The domain of integration of the n-fold integral is reduced to supp

ρ 1 × • • • × supp ρ n .
We first prove that the measure μ =

σ∈Σ n (σ)µ σ(1) ⊗ • • • ⊗ µ σ(n) is positive on the product space (supp ρ 1 , L)ו • •×(supp ρ n , L)
equipped with the product σ-field (where L denotes the one-dimensional Lebesgue σ-field). Notice that the product σ-field L ⊗n does not coincide in general with the n-dimensional Lebesgue σ-field (i.e. the completion of the n-dimensional Borel σ-field). Consider first the case of a subset of supp

ρ 1 × • • • × supp ρ n which is a product set A 1 × • • • × A n (where the A i 's are measurable). Necessarily, each A i is a subset of supp ρ i . If none of the A i 's is negligible, then we have A 1 < • • • < A n and μ(A 1 × • • • × A n ) = det[µ(A 1 ), • • • , µ(A n )]
is positive by definition. Suppose now some of the A i 's are negligible. For each index i, 1 ≤ i ≤ n, there exists a decreasing sequence (B i m ) m∈N of non-negligible measurable subsets of supp ρ i having an empty intersection (this is a consequence of the fact that supp ρ i is not negligible). Now for each m we have

A 1 ∪ B 1 m < • • • < A n ∪ B n m whence μ(A 1 ∪ B 1 m × • • • × A n ∪ B n m ) is positive. By the continuity of the measure µ we have μ(A 1 × • • • × A n ) = lim m→∞ μ(A 1 ∪ B 1 m × • • • × A n ∪ B n m ) so that μ(A 1 • • • A n ) is non-negative.
It follows that μ is non-negative on the boolean algebra of the finite (disjoint) union of product sets: its unique extension to the σ-field L ⊗n generated by these products is also non-negative.

The function (s

1 , • • • , s n ) → ρ 1 (s 1 ) • • • ρ n (s n
) is positive everywhere on this set and is measurable with respect to the σ-field L ⊗n : thus the integral ρ 1 (s 1 )

• • • ρ n (s n ) dμ(s 1 , • • • , s n ) is positive.
Remark. If µ is absolutely continuous with respect to the Lebesgue measure then Lyapunov theorem yields an alternative proof of theorem 2.2. In fact

∀k ∈ {1, • • • , n} ∃E k ⊂ supp ρ k µ(ρ k ) = µ(E k ).
Necessarily µ(E k ) is non-zero for each k (see remark after definition 2.1) and the absolute continuity hypothesis on µ implies that the E k 's are not negligible.

It follows that E 1 < • • • < E n and det[µ(ρ 1 ), • • • , µ(ρ n )] = det[µ(E 1 ), • • • , µ(E n )] > 0.
We shall denote by Γ k the subset

Γ k = { (x 1 , • • • , x k ) ∈ R k : a ≤ x 1 ≤ • • • ≤ x k ≤ b }.
Definition 2.4. The measure µ is said to be locally oriented if for each n-tuple x of Γ n there exists a neighbourhood

V = V 1 × • • • × V n of x such that for each k-tuple of measurable sets A 1 < • • • < A k satisfying A 1 × • • • × A k ⊂ V 1 × • • • × V k , the determinant µ 1 (A 1 ) • • • µ 1 (A k ) . . . . . . . . . µ k (A 1 ) • • • µ k (A k ) is positive.
As a curiosity, we prove the following Proof. Let µ be a locally oriented measure. The compact set Γ n can be covered by a finite family of open sets (V i ) i∈Υ where

V i = I i 1 × • • • × I i n and (I i k ) i∈Υ 1≤k≤n are subintervals of [a, b] in such a way that for each k-tuple of measurable sets A 1 < • • • < A k satisfying A 1 × • • • × A k ⊂ V i for some i ∈ Υ, the determinant formed with the first k components of the vectors µ(A 1 ), • • • , µ(A k ) is positive.
Let (J l ) l∈Σ be the finite family of the atoms of the algebra generated by the sets (

I i k , i ∈ Υ, 1 ≤ k ≤ n) (thus the J l 's are exactly the sets of the form i,k:x∈I i k I i k for some x ∈ [a, b]). Let us remark that for each (l 1 , • • • , l k ) in Σ k , the product J l 1 × • • • × J l k is contained in some product I i 0 1 × • • • × I i 0 k . In fact J l 1 × • • • × J l k ⊂ i∈Υ I i 1 × • • • × I i k so that there exits i 0 such that J l 1 × • • • × J l k ∩ I i 0 1 × • • • × I i 0 k is not empty. It follows that J l 1 ∩ I i 0 1 = ∅,• • • , J l k ∩ I i 0 k = ∅
and by the very construction of the sets J l 's we obtain

J l 1 ⊂ I i 0 1 , • • • , J l k ⊂ I i 0 k . We denote by μk the measure μk = σ∈Σ k (σ)µ σ(1) ⊗ • • • ⊗ µ σ(k) . Let (A 1 , • • • , A k ) be a k-tuple of measurable sets such that A 1 < • • • < A k . The product A 1 ו • •×A k is the disjoint union of the sets (A 1 ו • •×A k )∩(J l 1 ו • •×J l k ) when (l 1 , • • • , l k ) varies in Σ k . Let now (l 1 , • • • , l k ) belong to Σ k . Either (A 1 × • • • × A k ) ∩ (J l 1 × • • • × J l k
) is empty (and thus has a zero μk measure) or it is not empty and necessarily,

J l 1 < • • • < J l k .
Proceeding as in the proof of theorem 2.2, we show that μk is a positive measure on the set (J

l 1 × • • • × J l k ) whence μk ((A 1 × • • • × A k ) ∩ (J l 1 × • • • × J l k )) is non-negative. Since the set A 1 × • • • × A k is not negligible, at least one of these sets is not negligible. Let (A 1 × • • • × A k ) ∩ (J l 1 × • • • × J l k ) be such a set. It's a subset of one of the V i 's and moreover (A 1 ∩ J l 1 ) < • • • < (A k ∩ J l k ) whence μk ((A 1 ∩ J l 1 ) × • • • × (A k ∩ J l k )) is positive. Thus μk (A 1 × • • • × A k ) is positive.

Oriented measures with densities

Orientation condition ∆. We say that n real functions

f 1 , • • • , f n verify condition ∆ on an interval [a, b] if for each k in {1, • • • , n}, the determinant f 1 (x 1 ) • • • f 1 (x k ) f 2 (x 1 ) • • • f 2 (x k ) . . . . . . . . . f k (x 1 ) • • • f k (x k ) is positive whenever the x i 's in [a, b] are such that a ≤ x 1 < x 2 < • • • < x k ≤ b.
Remark. In our previous paper [START_REF] Cerf | Oriented measures with continuous densities and the bang-bang principle[END_REF], we didn't impose the sign of the above determinant to be positive. When dealing with continuous functions, a connectedness argument shows immediately that the sign is constant on the set Γ k . In our present framework (at the measure level), we find it convenient to work with this slightly more restrictive condition.

Examples. For n = 1, condition ∆ states that the function f 1 is positive; for n = 2, the functions f 1 , f 2 satisfy ∆ if and only if f 1 is positive and f 2 /f 1 is strictly increasing. The functions f i (t) = t i-1 (i ≥ 1) satisfy condition ∆ on R (the corresponding determinants are Vandermonde determinants).

= (µ 1 , • • • , µ n ) is oriented. Proof. Let A 1 < • • • < A k be k measurable sets of [a, b].
Since the determinant is a multilinear continuous form, we can write

A 1 f 1 • • • A k f 1 . . . . . . . . . A 1 f k • • • A k f k = • • • A 1 ו••×A k f 1 (s 1 ) • • • f 1 (s k ) f 2 (s 1 ) • • • f 2 (s k ) . . . . . . . . . f k (s 1 ) • • • f k (s k ) ds 1 • • • ds k .
By condition ∆, the integrand is positive on

A 1 × • • • × A k . If f 1 , • • • , f k are of class C k-1 on [a, b] we will denote their Wronskian by W (f 1 , • • • , f k ).
The following operational criterion for the fulfilment of the orientation condition ∆ has been used in [START_REF] Cerf | Oriented measures with continuous densities and the bang-bang principle[END_REF].

Proposition 3.2. Let f 1 , • • • , f n ∈ C n-1 ([a, b]) be such that ∀t ∈ [a, b] W (f 1 )(t) > 0, • • • , W (f 1 , • • • , f n )(t) > 0.
Then

f 1 , • • • , f n satisfy the orientation condition ∆ on [a, b].

Notations and preliminary lemmas

Let us introduce some notations.

If u 1 , • • • , u n are vectors of R n , their determinant is sometimes denoted by det [u 1 , • • • , u n ].
Let A be a n × n matrix with real coefficients; by det A or |A| we denote its determinant.

For each i, j ∈ {1, • • • , n}, by A ij we mean the (n -1) × (n -1) matrix obtained by removing the i-th row and the j-th column from A. Surprisingly, the following simple algebraic trick will play an essential role in the existence part of the proof of theorem 1.

Lemma 4.1. Let A = (a ij ) 1≤i,j≤n be an n×n matrix with real coefficients. Let

x 1 , • • • , x n be such that        a 1,1 x 1 + • • • + a 1,n-1 x n-1 + a 1,n x n = 0 a 2,1 x 1 + • • • + a 2,n-1 x n-1 + a 2,n x n = 0 . . . . . . . . . . . . a n-1,1 x 1 + • • • +a n-1,n-1 x n-1 +a n-1,n x n = 0 If det A nn = 0 then a n1 x 1 + • • • + a nn x n = |A| |A nn | x n .
Proof. Cramer rule applied to the above system yields

∀i ∈ {1, • • • , n -1} x i = (-1) n+i |A ni | |A nn | x n so that a n1 x 1 + • • • + a nn x n = n i=1 (-1) n+i a ni |A ni | |A nn | x n = |A| |A nn | x n since n i=1 (-1) n+i a ni |A ni | is the development of the determinant |A| along the first row.
The next lemmas involve strongly the notion of oriented measure.

F = l i=1 I i , G = m j=1 J j and let µ = (µ 1 , • • • , µ n ) be an oriented measure. Assume µ(F ) = µ(G). Then n < l + m; moreover if ∂F ∩ ∂G = ∅ then n < l + m -1.
Proof. Let us first remark that the symmetric difference

(I 1 ∪ • • • ∪ I l ) ∆ (J 1 ∪ • • • ∪ J m ) = i,j (I i ∪ J j ) \ i,j (I i ∩ J j )
is the union of at most l + m non-trivial intervals and that whenever at least two intervals have a common boundary point then this number is smaller than l + m -1. Since the intervals I 1 , • • • , I l are disjoint, as well as J 1 , • • • , J m , we have

(I 1 ∪• • •∪I l )∪(J 1 ∪• • •∪J m )\(I 1 ∩J 1 ) = (I 1 ∪J 1 )\(I 1 ∩J 1 )∪(I 2 ∪• • •∪I l )∪(J 2 ∪• • •∪J m ). Now, the set (I 2 ∪ • • • ∪ I l ) ∪ (J 2 ∪ • • • ∪ J m ) is a union of at most l + m -2 disjoint intervals.
Either I 1 ∩J 1 = ∅ or I 1 ∩J 1 = ∅ and (I 1 ∪J 1 ) is an interval. In both cases (I 1 ∪J 1 )\(I 1 ∩J 1 ) is the union of at most two intervals (at most one if I 1 and J 1 have a boundary point in common). A straightforward induction gives the result. Since the sets F and G are distinct, F ∆G is not empty. Let A 1 < • • • < A p be the connected components of F ∆G. For k in {1, • • • , p} we have

A k = (A k ∩ F ) ∪ (A k ∩ G) , (A k ∩ F ) ∩ (A k ∩ G) ⊂ A k ∩ (F ∩ G) ⊂ (F ∆G) ∩ (F ∩ G) = ∅ ; the set A k being connected, either A k ⊂ F \ G or A k ⊂ G \ F . Put λ k = +1 if A k ⊂ F \ G -1 if A k ⊂ G \ F
so that the equality µ(F ) = µ(G) can be rewritten as

     λ 1 µ 1 (A 1 )+ • • • +λ p µ 1 (A p )= 0 . . . . . . . . . λ 1 µ n (A 1 )+ • • • +λ p µ n (A p )= 0
Suppose n ≥ p; the first p equations imply that the determinant

µ 1 (A 1 ) • • • µ 1 (A p ) . . . . . . . . . µ p (A 1 ) • • • µ p (A p )
vanishes, which contradicts the fact that µ is oriented.

The following notations will be used throughout the paper.

Notations 4.3. We shall denote by Γ k the set

Γ k = { (γ 1 , • • • , γ k ) ∈ R k : a ≤ γ 1 ≤ • • • ≤ γ k ≤ b }.
To each k-tuple γ = (γ 1 , • • • , γ k ) belonging to Γ k we associate the two sets

E - γ = 0≤i≤k i odd [γ i , γ i+1 ] , E + γ = 0≤i≤k i even [γ i , γ i+1 ]
where by convention γ 0 = a, γ k+1 = b.

= (γ 1 , • • • , γ n ) and δ = (δ 1 , • • • , δ n ) of Γ n satisfy µ(E - γ ) = µ(E - δ ) (respectively µ(E + γ ) = µ(E + δ )). Then E - γ = E - δ (resp. E + γ = E + δ ).
Proof. Assume E - γ , E - δ are distinct and µ(E - γ ) = µ(E - δ ). Now, two possible cases may occur according to the parity of n.

• If n = 2r the sets E - γ and E - δ are the union of at most r intervals. Lemma 4.2 implies n < r + r which is absurd.

• If n = 2r + 1 the sets E - γ and E - δ are the union of at most r + 1 intervals. However b is a common boundary point. Lemma 4.2 implies n < (r + 1) + (r + 1) -1 which is absurd. The dual case µ(E + γ ) = µ(E + δ ) can be treated similarly. The following essential lemma will be used repeatedly. Proof. Consider the n × n linear system

λ 0 µ(I 0 ) -λ 1 µ(I 1 ) + • • • + (-1) n-1 λ n-1 µ(I n-1 ) = (-1) n-1 λ n µ(I n ).
where λ n is a parameter. The determinant of the system is

ω n = (-1) n(n-1) 2 det [µ(I 0 ), • • • , µ(I n-1 )].
The measure µ being oriented, ω n is not zero. Moreover, for each i in {0, • • • , n -1},

ω i = µ 1 (I 0 ) • • • (-1) i-2 µ 1 (I i-2 ) (-1) n-1 µ 1 (I n ) (-1) i µ 1 (I i ) • • • (-1) n-1 µ 1 (I n-1 ) µ 2 (I 0 ) • • • (-1) i-2 µ 2 (I i-2 ) (-1) n-1 µ 2 (I n ) (-1) i µ 2 (I i ) • • • (-1) n-1 µ 2 (I n-1 ) . . . . . . . . . . . . . . . . . . . . . µ n (I 0 ) • • • (-1) i-2 µ n (I i-2 ) (-1) n-1 µ n (I n ) (-1) i µ n (I i ) • • • (-1) n-1 µ n (I n-1 )
i.e. ω i = (-1)

n(n-1) 2 det [µ(I 0 ), • • • , µ(I i-2 ), µ(I i ), • • • , µ(I n )].
By Cramer formula, λ i equals λ n ω i /ω n . The measure µ being oriented ω i and ω n have the same sign so that λ i is positive whenever λ n is positive. Choosing λ n such that

0 < λ n < min( ω n ω 0 , • • • , ω n ω n-1 , )
we obtain an (n + 1)-tuple which solves the problem.

Main result

The statement of the main result involves the notations 4.3.

= (α 1 , • • • , α n ) and β = (β 1 , • • • , β n ) in Γ n such that µ(E - α ) = b a ρ dµ = µ(E + β ). ( * )
If in addition 0 < ρ < 1 then α and β in Γ n satisfying ( * ) are unique and verify

a < α 1 < • • • < α n < b , a < β 1 < • • • < β n < b.
Remark. The measure µ being non-atomic we don't care about boundary points of intervals and we might write µ(α, β) for the measure of the interval µ([α, β]).

Proof. We consider first the case 0 < ρ < 1 and we prove the result by induction on n.

µ([α 1 , b]) = b a ρ dµ = µ([a, β 1 ]).
• Assume the result is true at rank n -1. We deal only with the n-tuple β: existence of the n-tuple α corresponding to ρ at rank n follows from the fact that it coincides with the n-tuple β corresponding to 1 -ρ.

Define for each k in {1, • • • , n} µ k (ρ) =
a < ᾱ1 < • • • < ᾱn-1 < b , a < β1 < • • • < βn-1 < b and for each k in {1, • • • , n -1} θ k (a, ᾱ1 , • • • , ᾱn-1 ) = 0≤i≤n-1 i odd µ k (ᾱ i , ᾱi+1 ) = µ k (ρ) , θ k ( β1 , • • • , βn-1 , b) = 0≤i≤n-1 i even µ k ( βi , βi+1 ) = µ k (ρ) . ( * * ) Put S = β = (β 1 , • • • , β n ) ∈ Γ n : β 1 ≤ β1 , ∀k ∈ {1, • • • , n -1} θ k (β) = µ k (ρ) . Since ( β1 , • • • , βn-1 , b
) and (a, ᾱ1 , • • • , ᾱn-1 ) belong to S, the set S is not empty.

We show now that

either θ n ( β1 , • • • , βn-1 , b) < µ n (ρ) < θ n (a, ᾱ1 , • • • , ᾱn-1 ) or θ n (a, ᾱ1 , • • • , ᾱn-1 ) < µ n (ρ) < θ n ( β1 , • • • , βn-1 , b).
The equalities ( * * ) yield for each

k in {1, • • • , n -1} 0≤i≤n-1 i even βi+1 βi (1 -ρ) dµ k - 0≤i≤n-1 i odd βi+1 βi ρ dµ k = 0. Put for k, j in {1, • • • , n} x β j = (-1) j+1 , a β kj = βj βj-1 ρ β j dµ k , A β = a β kj 1≤k,j≤n
where

ρ β j = ρ if j is even, 1 -ρ if j is odd.
With these notations the above equalities become

∀k ∈ {1, • • • , n -1} n j=1 a β kj x β j = 0.
Since the measure µ is oriented then the determinant |A β nn | does not vanish by theorem 2.2. We are thus in the position to apply lemma 4.1:

θ n ( β1 , • • • , βn-1 , b) -µ n (ρ) = n j=1 a β nj x β j = |A β | |A β nn | (-1) n+1 .
Similarly if we define for k, j in {1, • • • , n}

x α j = (-1) j , a α kj = ᾱj ᾱj-1 ρ α j dµ k , A α = a α kj 1≤k,j≤n
where

ρ α j = ρ if j is odd, 1 -ρ if j is even, we have θ n (a, ᾱ1 , • • • , ᾱn-1 ) -µ n (ρ) = |A α | |A α nn | (-1) n .
The

, • • • , β n ) satisfying for each k in {1, • • • , n -1} µ k (a, β 1 ) + 2≤i≤n i even µ k (β i , β i+1 ) = µ k (ρ) = µ k (a, β1 ) + 2≤i≤n-1 i even µ k ( βi , βi+1 ) or equivalently ∀k ∈ {1, • • • , n -1} 2≤i≤n i even µ k (β i , β i+1 ) = µ k (β 1 , β1 ) + 2≤i≤n-1 i even µ k ( βi , βi+1 )
Suppose first β 1 = β1 . The above equations become

∀k ∈ {1, • • • , n -1} 2≤i≤n i even µ k (β i , β i+1 ) = 2≤i≤n-1 i even µ k ( βi , βi+1 ). We put β = (β 2 , • • • , β n-1 , β n ) and β = ( β2 , • • • , βn-1 , b).
If n is odd then

E - β = [β 2 , β 3 ] ∪ • • • ∪ [β n-1 , β n ], E - β = [ β2 , β3 ] ∪ • • • ∪ [ βn-1 , b];
if n is even then

E - β = [β 2 , β 3 ] ∪ • • • ∪ [β n , b], E - β = [ β2 , β3 ] ∪ • • • ∪ [ βn-2 , βn-1 ].
In both cases the preceding formulae can be rewritten as

∀k ∈ {1, • • • , n -1} µ k (E - β ) = µ k (E - β ); lemma 4.4 implies that E - β = E - β . Since in addition β2 < • • • < βn-1 < b then necessarily β 2 = β2 , • • • , β n-1 = βn-1 , β n = b. Suppose now β < β1 . Since β 1 < β1 < • • • < βn-1 < b then lemma 4.5 yields the existence of n real numbers λ 1 , • • • , λ n in ]0, 1/2[ such that for each k in {1, • • • , n -1} -λ 1 µ k (β 1 , β1 ) + 1≤i≤n-1 (-1) i+1 λ i+1 µ k ( βi , βi+1 ) = 0. The function ρ = (1 -λ 1 )χ [β 1 , β1 ] + 1≤i≤n-1 i odd λ i+1 χ [ βi , βi+1 ] + 2≤i≤n-1 i even (1 -λ i+1 )χ [ βi , βi+1 ] satisfies 0 < ρ < 1 on [β 1 , b] and for each k in {1, • • • , n -1} b β 1 ρ dµ k = µ k (β 1 , β1 ) + 2≤i≤n-1 i even µ k ( βi , βi+1 ).
We are thus led to find a (n -1)-tuple (

β 2 , • • • , β n ) such that (β 1 ≤)β 2 ≤ • • • ≤ β n (≤ b) and for each k in {1, • • • , n -1} 2≤i≤n i even µ k (β i , β i+1 ) = b β 1 ρ dµ k , or equivalently, if we put β = (β 2 , • • • , β n ), ∀k ∈ {1, • • • , n -1} µ k (E - β ) = b β 1 ρ dµ k .
Existence and uniqueness of β follow from the inductive assumption at rank n -1.

In addition, since 0

< ρ < 1 on [β 1 , b], we have β 1 < β 2 < • • • < β n < b.
We can thus define a map ψ :

[a, β1 ] → R n-1 such that for all n-tuple (β 1 , • • • , β n ) in Γ n (β 1 , • • • , β n ) ∈ S ⇐⇒ (β 2 , • • • , β n ) = ψ(β 1 ).
Thus S is the graph of ψ. By the continuity of the measure µ, the maps θ k , 1 ≤ k ≤ n -1, are continuous so that the set S is closed; moreover, the function ψ takes its values in the compact set [a, b] n-1 . It follows that ψ is continuous. Henceforth S is connected. As a consequence, the map θ n , being continuous on S, reaches all the values between θ n ( β1 , • • • , βn-1 , b) and θ n (a, ᾱ1 , • • • , ᾱn-1 ). In particular, there exists a n-tuple β in S such that θ n (β) = µ n (ρ). This n-tuple β solves the problem. Since θ n (a, ᾱ1 ,

• • • , ᾱn-1 ) = µ n (ρ) and θ n ( β1 , • • • , βn-1 , b) = µ n (ρ) then a < β 1 < β1 so that a < β 1 < β 2 < • • • < β n < b.
Uniqueness of β follows from lemma 4.4.

Consider now the case 0 ≤ ρ ≤ 1. Let (ρ m ) m∈N be a sequence of measurable functions such that 0 < ρ m < 1 and ρ m converges to ρ in L 1 µ ([a, b]). For each function ρ m there exists a unique n-tuple β m such that

µ(E + β m ) = b a ρ m dµ.
By compactness, we may assume that β m converges to some n-tuple β of Γ n . Passing to the limit, we obtain µ(E + β ) = µ(ρ). Finally we show that the map θ is proper (i.e. that the inverse image of a compact subset is compact). Let K be a compact subset of F and (γ m ) m∈N be a sequence in θ -1 (K) such that θ(γ m ) converges to µ(ρ) for some ρ, 0 < ρ < 1. Since the sequence (γ m ) m∈N is contained in Γ n , by compactness, we may assume that γ m converges to γ in Γ n . By the continuity of θ, we have

θ(γ) = µ(E - γ ) = b a ρ dµ.
Uniqueness part of theorem 5.1 implies that γ belongs to Γ • n . The map θ is proper and thus closed. It follows that its inverse θ -1 is continuous. The equality θ(∂Γ n ) = ∂R is a consequence of the inclusion θ(Γ

• n ) ⊂ R
• and the fact that θ is one to one.

We refer to [START_REF] Valentine | Convex sets[END_REF] for the definitions of classical notions associated with convex sets. We have the following Theorem 6.5. The range R of an oriented measure is strictly convex.

Proof. Let µ(E), µ(F ) be two distinct points of R. By theorem 5.1 we may assume that the sets E and F are finite unions of closed intervals. Let λ ∈]0, 1[ and put ρ = λχ E +(1-λ)χ F . Assume for instance E \ F = ∅. Then there exists a non-trivial interval I such that ∀x ∈ I ρ(x) = λχ E (x) + (1 -λ)χ F (x) = λ.

Put = min(λ, 1 -λ). Lemma 6.1 applied to ρ, I, shows that µ(ρ) belongs to R • .

  Lyapunov states that the range R of a non-atomic vector measure µ on [a, b] R = { µ(A) : A measurable subset of [a, b] } coincides with the convex set b a ρ dµ : 0 ≤ ρ ≤ 1 .

Proposition 2 . 5 .

 25 A locally oriented measure on [a, b] is oriented on [a, b].

Proposition 3 . 1 .

 31 Let f 1 , • • • , f n be n functions in L 1 ([a, b]) satisfying the orientation condition ∆ on [a, b]. Let µ i be the measure on [a, b] whose density with respect to the Lebesgue measure is f i . Then the measure µ

Lemma 4 . 2 .

 42 Let F and G be two distinct subsets of [a, b] which are the union of l and m disjoint closed intervals

Lemma 4 . 4 (

 44 Uniqueness). Let µ be a n-dimensional oriented measure on [a, b]. Assume the n-tuples γ

Lemma 4 . 5 .(- 1 )

 451 Let µ = (µ 1 , • • • , µ n ) be an oriented measure on the interval [a, b] and I 0 < I 1 < • • • < I n be n + 1 subintervals of [a, b]. Then, given a positive , there exist n + 1 positive real numbers λ 0 , • • • , λ n such that ∀l ∈ {0, • • • , n} 0 < λ l < and n l=0 l λ l µ(I l ) = 0.

Theorem 5 . 1 .

 51 Let µ be an oriented measure on [a, b] and let ρ be a measurable function defined on [a, b] with values in [0, 1]. There exist two n-tuples α

  • n=1. The measure µ being oriented on [a, b], the maps α → µ([α, b]) and β → µ([a, β]) are continuous and strictly monotonic on [a, b]. It follows that there exist unique real numbers α 1 and β 1 such that

  and for each n-tuple β in Γ n θ k (β) = µ k (E + β ). The inductive assumption yields the existence of two (n -1)-tuples ᾱ = (ᾱ 1 , • • • , ᾱn-1 ) and β = ( β1 , • • • , βn-1 ) such that

6 .

 6 The range of an oriented measure Let µ be an oriented measure on [a, b]. We denote by R the range of µ i.e. R = { µ(A) : A measurable subset of [a, b] }. Lemma 6.1. Let ρ be a measurable function on [a, b], 0 ≤ ρ ≤ 1. Suppose there exist a non-trivial interval I of [a, b] and a positive real number such that ≤ ρ ≤ 1 -on I. Then the set b a ρ dµ : ρ= νχ I + ρ, ν ∈ L 1 µ (I), |ν| < is a neighbourhood of b a ρ dµ in R n . Proof. Let I 1 < • • • < I n be n non-trivial subintervals of I. The measure µ being oriented, the vectors µ(I 1 ), • • • , µ(I n ) form a basis of R n . The map Λ : (λ 1 , • • • , λ n ) ∈ R n -→ 1≤i≤n λ i µ(I i ) ∈ R n is a linear isomorphism and is thus open. Let V = { (λ 1 , • • • , λ n ) : max 1≤i≤n |λ i | < }.Since Λ(V ) is a neighbourhood of the origin and is contained in the setI ν dµ : ν ∈ L 1 µ (I), |ν| < ,piecewise constant function ρ such that 0 < ρ < 1 and µ(ρ) = θ(γ). Clearly there exist a positive and a subinterval I of [a, b] on which ≤ ρ ≤ 1 -. PutV I, ρ = { νχ I + ρ : ν ∈ L 1 µ (I), |ν| < }.Lemma 6.1 implies that the setµ(V I, ρ ) = b a ρ dµ : ρ ∈ V I, ρ is a neighbourhood of µ(ρ) in R n . Since each element ρ of V I, ρ satisfies 0 < ρ < 1 then µ(V I, ρ ) is entirely contained in F . Moreover F coincides with θ(Γ • n ) and thus θ(Γ • n ) is a neighbourhood of θ(γ).Now each open convex set in R n is the interior of its closure; by lemma 6.3, the set θ(Γ • n ) is convex and its closure is R, whence θ(Γ• n ) = R • .

  measure µ being oriented, the determinants |A α | and |A β | have the same sign, as do |A α nn | and |A β nn |. It follows that θ n ( β1 , • • • , βn-1 , b)-µ n (ρ) and θ n (a, ᾱ1 , • • • , ᾱn-1 )-µ n (ρ) have opposite signs.At this stage, we prove that the set S is the graph of a continuous function, this will imply that S is connected. Let β 1 belong to [a, β1 ]. We are looking for a (n -1)-tuple (β 2

thank the D.M.I. of the E.N.S. for the technical support during the preparation of this paper. The second author (C.M.) was partially supported by a Grant of the Consiglio Nazionale delle Ricerche (Grant 203.01.62).

then the conclusion follows.

Remark. The hypothesis ≤ ρ ≤ 1 -implies that µ(ρ) belongs to the interior of R.

Remark. The conclusion of lemma 6.1 does not hold for an arbitrary vector measure: consider for instance the n-dimensional Lebesgue measure.

Let θ : Γ n → R be the function defined by θ(γ) = µ(E - γ ). The interior of Γ n is the set Γ

• n ; applying lemma 4.5 to µ, γ and < 1/2, we obtain a (n + 1)-tuple (λ 0 ,

.

By construction we have 0 < ρ < 1 and

so that θ(γ) belongs to F .

We have the following Theorem 6.4. The range of θ coincides with R; the map θ induces an homeomorphism from Γ • n onto R • and maps ∂Γ n onto ∂R.

Proof. The surjectivity of θ follows directly from theorem 5.1. Injectivity of the restriction of θ to Γ • n is a consequence of the uniqueness part of theorem 5.1 together with lemma 6.3. We claim that θ(Γ Proof. We first remark that the family of the sets which are a finite union of intervals and whose characteristic function has less than n -1 discontinuity points coincides with the family { E - γ : γ ∈ ∂Γ n }. Theorem 6.4 shows that µ(F ) belongs to ∂R whenever F = E - γ for some γ ∈ ∂Γ n . Conversely let E be such that µ(E) belongs to ∂R. Theorem 6.4 yields the existence of a n-tuple γ belonging to ∂Γ n such that µ(E - γ ) = µ(E); a consequence of theorem 6.5 is that µ(E) is an extreme point of R. Olech Theorem [5, Th. 1] implies that E∆E - γ is µ-negligible.

Our approach discloses the recursive structure of the boundary of the range of an oriented measure. For k belonging to {0, • • • , n} let

Proof. Injectivity follows directly from corollary 6.6. The rest of the proof uses the techniques of the proof of theorem 6.4.

These results yield the following Proposition 6.8. The boundary of the range R of an oriented n-dimensional measure admits the decomposition

Let T be the symmetry with respect to µ(a, b)/2 (so that for each measurable subset A of [a, b], T (µ(A)) = µ([a, b] \ A)). Then for each k belonging to {0, • • • , n} we have


S Barthélémy

M Parisot

M.-H Le D,

N Goutal

S Ricci

Lateral coupling method for different one and two dimensional shallow-water solvers: application to the hydro-informatic Telemac-Mascaret system

Keywords: shallow-water equations, lateral coupling, inundation simulation, river channel, floodplain, finite volume, finite element, Telemac-Mascaret, OpenPALM

We present a lateral coupling method between one-dimensional (1D) and two-dimensional (2D) shallow-water solvers dedicated to numerical simulation of river overflows. Both 1D and 2D models, respectively in the river channel and in the floodplains, can be solved with implicit scheme in order to limit computational cost. The 1D-2D exchange terms for mass and momentum between the river and the floodplains are computed with an explicit solver of Riemann problem over a local area that extends over the coupling interface. A distinctive feature of the approach is its flexibility in reusing available computation codes of 1D and 2D shallow-water models as black-boxes. The proposed method has been implemented with the integrated suite of solvers Telemac-Mascaret together with the dynamical coupling software OpenPALM. Numerical validations with respect to an analytical solution and experimental data as well as a first application for the historical 1981 flood event over the Garonne River are presented. The results showed that the overflow discharge is well estimated by the coupled model while it tends to underestimate the water height near the lateral boundaries. For the real application case, reasonable agreement was found between the full 2D simulation and the coupled model.

Introduction

Shallow-water models are widely used in the field of rivers and maritime hydraulics. These equations are derived from the Navier-Stokes equations for shallow-flows assuming hydrostatic pressure and low variation in bathymetry. In the context of hydrodynamic modeling, the use of full twodimensional (2D) shallow-water equations is well adapted in areas where the behavior of flows is no longer mono-dimensional, such as in confluence zones or floodplains. Nevertheless extensive use of full 2D model for operational flood forecasting may be limited due to the lack of bathymetry data, numerical difficulties in the case of steep slopes (e.g. by dyke system) and also due to computational cost constraints.

Multi-dimensional coupling offers an appropriate solution with one-dimensional (1D) model where the flow is nearly uni-directional and with local 2D models where needed. This solution makes the most of the benefits of 1D models and allows for relevant representation of complex processes with 2D models. The computational cost of such coupling solution is significantly smaller than that of a full 2D simulation, while maintaining better simulated hydraulic state than that of the full 1D model, see. e.g. [START_REF] Barthélémy | On operational flood forecasting system involving 1d/2d coupled hydraulic model and data assimilation[END_REF][START_REF] Morales-Hernández | Conservative 1d-2d coupled numerical strategies applied to river flooding: The tiber (rome)[END_REF][START_REF] Malleron | On the use of a high-performance framework for efficient model coupling in hydroinformatics[END_REF].

When the river overflows into its surrounding plains or, inversely, when the plains empty into the river, lateral exchanges between the river and its floodplains occur. This situation can be modeled either by a full 2D model or by a lateral coupling between 1D section-averaged shallowwater equations in the river channel and 2D equations for the flow in the floodplain. This last strategy allows to keep the advantages of each model in its respective sub-domain, and to optimize the computation cost. Furthermore, it avoids the need of local 2D mesh refinement where strong gradients in bathymetry occur; for instance dikes can be included within the 1D sub-domain.

Different coupling strategies adapted to overflowing simulations are reported in the literature and can be classified into two main categories: overlapping and non-overlapping coupling. Pioneer works by [START_REF] Fernández-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF][START_REF] Gejadze | On a 2d 'zoom' for the 1d shallow water model: Coupling and data assimilation[END_REF][START_REF] Marin | Superposition of local zoom models and simultaneous calibration for 1d-2d shallow water flows[END_REF] present an overlapping strategy for which the 2D model is also solved in the river channel. Such a coupling strategy remains subject to previously listed difficulties (dry zones, strong gradient in bathymetry). In the non-overlapping coupling approach, 1D and 2D models are solved separately over their respective sub-domain; the coupling is achieved via a prescription of boundary fluxes for the 2D model, while injecting exchange terms of mass and momentum into the 1D model as source terms. The estimation of these coupling quantities is based on the resolution of a 2D Riemann problem at the lateral interfaces. This later coupling strategy succeeds in preserving some required characteristics of a coupled solvers such as conservativity, positivity and well-balancing; it is thus widely used [START_REF] Bladé | Integration of 1d and 2d finite volume schemes for computations of water flow in natural channels[END_REF][START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF][START_REF] Morales-Hernández | A conservative strategy to couple 1d and 2d models for shallow water flow simulation[END_REF][START_REF] Morales-Hernández | A riemann coupled edge (rce) 1d-2d finite volume inundation and solute transport model[END_REF][START_REF] Morales-Hernández | Conservative 1d-2d coupled numerical strategies applied to river flooding: The tiber (rome)[END_REF][START_REF] Bouchut | Nonlinear Stability of Finite Volume Methods for Hyperbolic Conservation Laws[END_REF].

The proposed coupling algorithm follows non-overlapping approach and presents two major advantages. First, it is non-intrusive to both 1D and 2D solvers. In the present work, the open source integrated suite of solvers Telemac-Mascaret1 is used to solve the shallow-water equations in 1D and 2D sub-domains; the 1D solver is based on a finite volume scheme and the 2D solver is based on a linear finite element scheme [START_REF] Hervouet | a hydroinformatic system[END_REF]. Second, it is compatible with large-timestep often used by implicit scheme of 2D solver. It means that no additional constraint on the Courant-Friedrichs-Lewy (CFL) criteria is required by the coupling. Following the work from [START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF], an explicit solver of 2D local Riemann problem is used to compute the exchange terms between 1D and 2D sub-domains. This procedure is confined to a local area covering the coupling interface -further called buffer zone.

The paper is organized as follows: Section 2 briefly recalls the 1D and 2D shallow-water equations in order to formulate of the lateral coupling source terms in the framework of finite volume or finite element schemes. Section 3 details the construction of the buffer zone and how it has been used in the coupling algorithm. Implementation of the lateral coupling strategy for an academic test case, is presented in Section 4. Validation and application of the proposed coupling to the historical 1981 flood event over the Garonne River are presented in Section 5. Conclusions and perspectives for further improvements of the method are given in Section 6.

Overview on lateral coupling for shallow-water models

Let's consider a free surface shallow-flows on a bathymetry featuring a river channel and its surrounding floodplains, as represented in Fig. 1a. The river is delimited by left and right banks represented with regular curves l(ξ) and r(ξ) where ξ is the hydraulic axis. More precisely, there exists a curvilinear coordinates couple (ξ, ζ)(x, y) such that the center of the river corresponds to the line ζ = 0, the left bank to the line ζ = l(ξ) and the right bank to the line ζ = r(ξ). Denoting Ω 1D and Ω 2D the sub-domains corresponding to the river (1D) and the floodplains (2D), we define Γ l,r the outward normal vectors to the left and right river banks of Ω 1D : The 1D and 2D shallow-water models together with the principles of non-overlapping lateral coupling are presented in the following.

Γ l ξ Ω 1D Ω 2D Γ r Right Floodplain Left Floodplain Ω 2D ζ (a) z 2D b z 1D b η 1D V 1D b V 2D b ζ b(ξ) (b)
Γ l := J ξζ ∂l ∂ξ

1D and 2D governing equations

Shallow-water equations are commonly used to describe the governing of shallow-flows on bathymetry with small slopes. This system of equations can be derived by depth-integrating the Navier-Stokes equations under the assumption of hydrostatic pressure. The conservative form of the 2D numerical model reads:

∂V ∂t + ∂F x ∂x + ∂F y ∂y = S b + S f (2.1)
where V(t, x) is the conservative variable and F x (V) = (F i x ) i=1,3 , F y (V) = (F i y ) i=1,3 are the fluxes defined as :

V :=   h hu hv   , F x (V) :=   hu hu 2 + gh 2 2 huv   , F y (V) :=   hv huv hv 2 + gh 2 2   .
Here, h(t, x, y) stands for the water height, u(t, x, y) = (u, v) stands for the depth-averaged horizontal velocity of the flow and g is the gravity constant. The source terms S b and S f on the right-hand side of (2.1) describe the variation in bathymetry and the bed friction as

S b (V, z b ) :=    0 -gh ∂z b ∂x -gh ∂z b ∂y    , S f (V) :=    0 -gh n 2 |u|u h 4/3 -gh n 2 |u|v h 4/3    ,
where z b (x, y) represents the river bed elevation and n is the Manning friction coefficient.

Shallow-flow on floodplains is generally two-dimensional while the flow in the river channel is rather uni-directional. Consequently, 1D modelling is relevant to describe the dynamics of river flows. Furthermore, a 1D model requires less geometry data and offers remarkable gains on computational resources compared with a full 2D simulation. Integrating the 2D equations (2.1) along a cross-section, from left bank l(ξ) to right bank r(ξ), formulates the 1D shallow-water model along the 1D hydraulic axis of the river -the curvilinear-axis ξ following the main direction of the flow. More details on model derivation are found in [START_REF] Fernández-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF][START_REF] Gejadze | On a 2d 'zoom' for the 1d shallow water model: Coupling and data assimilation[END_REF][START_REF] Marin | Superposition of local zoom models and simultaneous calibration for 1d-2d shallow water flows[END_REF]. The 1D equations read

∂W ∂t + ∂G ∂ξ = B - b={l,r} O b , (2.2)
where W(t, ξ) := (S, Q) and G(ξ, W) are the state variables and the flux respectively. The wetted area S of cross-section (also called hydraulic section [START_REF] Goutal | A finite volume solver for 1d shallow-water equations applied to an actual river[END_REF]), the discharge Q in the river channel direction n ξ = ( ∂ξ ∂x , ∂ξ ∂y ) and the flux G are defined by

S := r(ξ) l(ξ) h 1D dζ, Q := r(ξ) l(ξ) h 1D u • n ξ dζ, G(ξ, W) := Q Q 2 S + P (ξ, S) with P (ξ, S) := g 2 r(ξ) l(ξ) (h 1D ) 2 dζ and h 1D (ξ, ζ, S) := max (0, η(ξ, S) -z b (ξ, ζ)) .
The bathymetry is projected onto the local (ξ, ζ) space as

z b (ξ, ζ) = z b J -1 ξζ ξ ζ
. The water level η in the river is assumed to be constant over a cross-section; it is hence computed solving the following inverse problem that involves the hydraulic section in the non-linear problem:

Find η(ξ, S) ∈ R such that r(ξ) l(ξ) max(0, η -z b ) dζ = S.
In practice, the water level η and the pressure P are tabulated during an initialization step in order to reduce the computational time.

The source term B reads

B(ξ, W) =   0 -gSJ -g r(ξ) l(ξ) h 1D ∂ z b ∂ξ dζ  
and represents both the friction and the variation in the cross-sections. The expression of the 2D friction law in the section-averaged framework leads to the source term

J(W) = n 2 Q 2 S 2 R 4/3 h
, where R h is the hydraulic radius, i.e. the ratio of hydraulic section and wetted perimeter of the cross-section. The remaining source terms O b ∈ R 2 of (2.2) stands for the lateral mass and momentum exchanges between river and floodplains due to overflows as long as the water heights at the banks are non-zero. In other words, it is the projection of 2D fluxes in the outward normal direction Γ b to the bank b, then oriented following the river. As stated in [START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF], O b reformulates as

O b = R ξ F(V b )Γ b 1,2 , b ∈ {l, r} (2.3)
where the notation (A) 1,2 ∈ R 2 stands for the two first components of the vector A ∈ R 3 . Here, V b is the 2D state variable evaluated at the river bank additionally, the rotation matrix R ξ towards the river channel direction n ξ and the flux F(V) are given by:

R ξ =    1 0 0 0 ∂ξ ∂x ∂ξ ∂y 0 -∂ξ ∂y ∂ξ ∂x    , F(V) =    F 1 x F 1 y F 2 x F 2 y F 3 x F 3 y    .

Formulation of the conservative lateral coupling

With regard to the 2D model, the lateral coupling relies on the computation of the boundary condition at the lateral 1D-2D interface, more precisely the 2D fluxes at the river banks. Solving the 2D Riemann problem along the 1D-2D interface is an adequate and natural way to do so, since this approach directly inherits the conservation property of the finite volume method.In this context, the overflow source term used in the 1D model is deduced from the 2D fluxes at the river banks.

The 2D Riemann problem depends on the 2D states (V 1D b , z b . Following [START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF], it reads

V 1D b (ξ, W) = R -1 ξ h 1D b   1 Q S 0   (2.4)
where h 1D b := h 1D (ξ, b(ξ), S) is the water height at the bank. This simple transformation neglects the local behaviors of flow near the bank the water level and the velocity at the bank are assumed to be close to that of the cross section and the velocity at the bank is assumed to be oriented along the hydraulic axis ξ of the river. More advanced transformation may be used, as proposed in [START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF] based on the successive resolution of Riemann problems at the bank.

Next, the following inhomogeneous Riemann problem is solved to identify the solution V b at the interface b(ξ):

         ∂V ∂t + ∂F x (V) ∂x + ∂F y (V) ∂y = S b (V, z b ), (V, z b )(0, x, y) = (V 1D b , z 1D b ) if |ζ(x, y)| < |b(ξ)|, (V 2D b , z 2D b ) otherwise. (2.5)
Note that in (2.5), the friction source term from the full 2D system (2.1) was dropped, since it does not contribute to the mass and momentum fluxes across the bank. Overflow source terms are then computed using equation (2.3).

Nevertheless, the direct resolution for (2.5) is known to be costly; it is more practical to replace the exact solution V b by an approximate one usually called simple solver. In turn, exchanged mass and momentum across the bank are estimated using an appropriate non-conservative numerical flux. Several well-known numerical fluxes for 2D shallow-water model are used for practical applications. We refer to [START_REF] Ata | A weighted average flux (waf) scheme applied to shallow water equations for real-life applications[END_REF] for a more detailed discussion. In the present work, the overflow source terms are estimated using the homogeneous HLLC flux F b , [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF], augmented by a hydrostatic source term S b , [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Accordingly, the overflow source term O b is estimated as:

O b = R ξ F b (V 1D b , V 2D b ; Γ b ) + S b (V 2D b , V 2D b ; Γ b ) 1,2 , (2.6)
where the states V 1D b and V 2D b are computed from the hydrostatic reconstruction, S b is derived by balancing the flux gradient with the bathymetry term for hydrostatic equilibrium (see [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] for more details)

V 1D b = (h, hu, hv) 1D b with h 1D b = max(0, h 1D b + z 1D b -max(z 1D b , z 2D b )), V 2D b = (h, hu, hv) 2D b with h 2D b = max(0, h 2D b + z 2D b -max(z 1D b , z 2D b )), S b (V 2D b , V 2D b ; Γ b ) = 0, - g 2 (h 2 -h 2 ) 2D b Γ b .
(2.7)

Finite volume discretization of the coupled 1D and 2D models

Estimating overflow source terms, which relies on solving Riemann problems, naturally leads to the use of finite volume solvers for both 1D and 2D models. This kind of coupling method has been widely adopted in recent works, see e.g. [START_REF] Bladé | Integration of 1d and 2d finite volume schemes for computations of water flow in natural channels[END_REF][START_REF] Fernández-Nieto | Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes[END_REF][START_REF] Morales-Hernández | A conservative strategy to couple 1d and 2d models for shallow water flow simulation[END_REF][START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF]. For clarification purposes, let us briefly recall the finite volume discretization of the coupled model.

C j Ω 2D Ω 1D l(ξ) Γ l x y 0 ξ i+1/2 ξ i+3/2 ξ i-3/2 ξ i-1/2 ξ Ω b C i
Figure 2: Sketch of river channel Ω 1D and its left floodplain Ω 2D discretized by 1D and 2D meshes respectively. The shadow sub-domain Ω b ⊂ Ω 2D stands for the buffer zone along the lateral boundary l(ξ) that will be used for our coupling method.

The domain decomposition into Ω 1D and Ω 2D is presented in Fig. 2. T 1D denotes the 1D mesh of the river channel, as a sequence of control volumes

I i = [ξ i-1/2 , ξ i+1/2 ] with ξ i+1/2 = ξ i +ξ i+1
2 for i ∈ T 1D and ξ i an increasing sequence of points along the curvilinear-axis. For each I i , a function p i (ζ), with l(ξ i ) ≤ ζ ≤ r(ξ i ), is used to represent the associated cross-section bathymetry of the river. The floodplain is discretized over a 2D mesh, T 2D , constituted of polygonal control volumes C j for j ∈ T 2D . The resolution of the 2D mesh is usually much finer than that of the 1D in order to solve more complex flow dynamics. As a consequence, a 1D cell I i is most likely adjacent to several 2D cells C j . Inversely, it could happen that a 2D cell might be adjacent to more than one cell of the 1D mesh because the meshing procedures on the sub-domains Ω 1D and Ω 2D are performed independently from each others.

Next, let us introduce a timestep, ∆t, assumed to be constant for simplicity purpose. Time discretization is defined by providing an increasing time sequence t n ∈ R + . Volume-averaged quantities of the solutions of 1D and 2D models at a time t n are respectively defined by

W n i = 1 |I i | I i W(t n , ξ) dξ, V n j = 1 |C j | C j V(t n , x, y) dx dy.
Finite volume schemes aim to compute the updated states W n+1 i and V n+1 j at time t n+1 = t n + ∆t from the given states W n i and V n j . Such a method is derived by integrating governing systems (2.1) and (2.2) on the corresponding control volumes within the time interval [t n , t n+1 ]. In order to preserve some relevant steady solutions of shallow-water equations (such a scheme is said wellbalanced), it would be preferred to perform the discretization of both flux and some well-chosen source terms at cell's interfaces [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][START_REF] Goutal | A finite volume solver for 1d shallow-water equations applied to an actual river[END_REF][START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF].

For the 1D model, the resulting finite volume scheme writes

W n+1 i = W n i - ∆t |I i | G L i+1/2 -G R i-1/2 - ∆t |I i | j∈L 2D b,i |b ij |O ij (2.8)
where

G R i-1/2 = G R (W n i-1 , W n i ; p i-1 , p i ) and G L i+1/2 = G L (W n i , W n i+1 ; p i , p i+1
) approximate the flux G(ξ, W) together with the contribution of the source term B(ξ, W) at the right-and left-sides of the cell's interfaces respectively, i.e. the quantities of mass and momentum entering and leaving the control volume I i during a timestep ∆t. See e.g. [START_REF] Goutal | A finite volume solver for 1d shallow-water equations applied to an actual river[END_REF] for such an example of well-balanced scheme. Lateral coupling between river and floodplains is accounted for by the last term in equation (2.8). For each 1D control volume I i , we denote by L 2D b,i ⊂ T 2D the set 2D cells adjacent to I i , |b ij | the length of corresponding lateral boundary. Finally, O ij stands for the lateral exchanges in mass and momentum fluxes computed by (2.6) after solving the Riemann problem at the bank.

Performing in similar way, one can write the resulting finite volume scheme for 2D model as

V n+1 j = V n j - ∆t |C j | i∈E j |e ji |F(V n j , V n i , z j , z i ; n ji ) + ∆t S f (V n j , V n+1 j ) + ∆t |C j | i∈L 1D b,j |b ij | F b (V 1D i , V n j ; Γ b ) + S b (V n j , V n j ; Γ b ) , (2.9)
where E j denotes the set of internal and all non-lateral boundary edges of C j ; the shorthand notation F(V n j , V n i , z j , z i ; n ji ) stands for a non conservative numerical flux (which includes bathymetry source terms) leaving the control volume C j through the face e ji = C j ∩ C i and in outward direc-

tion n ij . Friction term S f (V n j , V n+1 j
) is usually approximated by a (semi-)implicit discretization in order to avoid additional restriction on the timestep. It should be noted that the second right hand term in (2.9) expresses the boundary condition between the 2D cell C j and all neighboring 1D cells denoted L 1D b,j ⊂ T 1D . Let us notice that we find again the numerical flux F b and the source term S b of equation (2.6) which formulate exactly O ij in (2.8), and therefore allows to couple the 1D and 2D models.

Lateral coupling strategy including buffer zone

Although lateral coupling methods employing (2.6), (2.8) and (2.9) are known to be straightforward and accurate, they suffer from two main limitations:

• first, both 1D and 2D schemes are derived in the framework of explicit finite volume methods.

This limits the ability to use solvers with finite different or finite element solvers. Furthermore, implicit schemes are often used in engineering software for stability and performance purposes;

• second, the overall coupling procedure is restricted under a same timestep, ∆t, which could become very severe due to the CFL condition imposed by the finite volume scheme. That is

∆t v max ∆x 2D ≤ 1, (3.1)
where ∆x 2D stands for the typical resolution of the 2D mesh, v max expresses the largest wave speed of Riemann problems at cell faces in the floodplains. As a consequence, fast overbank flooding over complex bathymetry could lead to very small timestep, see e.g. [START_REF] Nicholas | Numerical simulation of overbank processes in topographically complex floodplain environments[END_REF].

In the following, a new coupling strategy will be presented that allows to overcome the aforementioned difficulties. The present method relies on the definition of a buffer zone -a local domain along the river banks -over which an explicit finite volume scheme is applied.

Buffer zone and exchange source terms

Along the lateral boundary b ∈ {l, r} between the river Ω 1D and a floodplain Ω 2D , we consider a sub-domain Ω b ⊂ Ω 2D called buffer zone. The key idea consists in only applying an explicit scheme, such that (2.9), over Ω b in order to accurately estimate the overflow source terms by (2.6) while maintaining the 1D and 2D models as black-box with their own numerical solvers. Numerical procedure implemented for the buffer zone is called the bank model. Even though the coupling strategy could be applied to any existing 1D and 2D solvers, it is here implemented with the Telemac-Mascaret suite of numerical codes, that imply the explicit finite volume scheme (2.8) for the 1D model and an implicit linear finite element scheme for the 2D model.

As such, the 2D state V n is defined by prescribing its values at the vertices of triangle mesh T 2D . Furthermore, these two solvers are well-balanced as they preserve at least the steady states of lake-at-rest, also called the C-property [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF].

On the buffer zone, a vertex-centered finite volume discretization is used in order to be compatible with the finite element solver of the 2D model. It means that Ω b is partitioned into dual cells {C j } j∈V b , where V b denotes the set of vertices of T 2D belonging to Ω b ; each dual cell C j is built by joining the gravity centers and the midpoint of edges of the triangles surrounding vertex j (see again Fig. 2).

Let us give more detail on the coupling procedure. First of all, one can notice that the resolution of the 1D mesh is in practice much coarser than that of the 2D mesh, typically a hundred meters in river channel compared with several meters for floodplain. As a consequence, relevant 1D timestep ∆t 1D is about one or two orders of magnitude compared to 2D timestep ∆t 2D . Furthermore, the coupling timestep ∆t (used for the bank model) might be even smaller than ∆t 2D due to CFL condition (3.1). Therefore, one can consider the following simplified setting:

∆t 1D = M ∆t 2D = M (K∆t) (3.2)
with some M, K ∈ N. It means that one timeloop of 1D scheme corresponds to M timeloops of 2D scheme and M × K sub-timeloops on the buffer zone. Assume that at time t n = n∆t 1D , the states W n and V n are known in Ω 1D and Ω 2D respectively. To carry out a cycle of the coupling procedure, i.e. to advance forward from t n to the next time t n+1 , the 2D simulation should perform M timeloops from starting value V n,0 = V n .

Consider now a time t n,m = t n + m∆t 2D , 0 ≤ m ≤ M -1, at which the 2D state V n,m is available. The coupling procedure is sketched in Fig. 3. It starts by performing K sub-timeloops, 0 ≤ k ≤ K -1, of explicit finite volume scheme (2.9) on the buffer zone in order to construct (h, q) n,m+1 j the water height and the discharge per width unit at the lateral boundary. Next, the 2D simulation computes the updated states V n,m+1 .

Once all M timeloops are achieved, in one hand one can obtain the updated 2D state V n+1 = V n,M and in the other hand we compute the updated 1D states W n+1 by performing 1D simulation provided the source term O 1D b constructed by the bank model during its M × K sub-timeloops. Due to the setting (3.2), the bank model must perform K sub-timeloops at the times t n,m,k = t n,m + k∆t, 0 ≤ k ≤ K -1, from t n,m to t n,m+1 starting with initial condition

W n W n+1 1D timeloop 1D -Mascaret W n i∈L 1D b,j V n,m i∈E ext j V n,m,k j V n,m,k+1 j sub-timeloop Bank -2D Riemann solver V n,m V n,m+1 2D timeloop 2D -Telemac V n,m,0 j∈V b O 1D b (h, q) n,m+1 j
V n,m,0 j = V n,m j , ∀j ∈ V b .
Since the main objective of these sub-timeloops is to estimate the exchange source terms, the friction term can be neglected when solving Riemann problems on the buffer zone. Therefore, finite volume scheme (2.9) used by these sub-timeloops becomes

V n,m,k+1 j = V n,m,k j - ∆t |C j | i∈E j |e ji |F(V n,m,k j , V n,m,k i , z j , z i ; n ji ) + ∆t |C j | i∈L 1D b,j |b ij |O n,m,k ij . (3.3)
Recall that on the buffer zone, one can choose to approximate the flux and the source terms with F b -a homogeneous HLLC flux augmented by the hydrostatic source term S b (2.7). Accordingly,

F(V n,m,k j , V n,m,k i , z j , z i ; n ji ) = F b (V n,m,k j , V n,m,k i ; n ji ) + S b (V n,m,k j , V n,m,k j ; n ji ), (3.4) O n,m,k ij = F b (V 1D (W n i ), V n,m,k j ; Γ b ) + S b (V n,m,k j , V n,m,k j ; Γ b ). (3.5)
Let us notice that the reconstructed states V n,m,k i for all i ∈ E ext j -the set of vertices of T 2D neighboring to C j and located outside of the buffer zone Ω b -remain not yet determined. Indeed, the 2D states V n,m i , for i ∈ E ext j , have not been updated during these sub-timeloops; the reconstructed states V n,m,k i do play the role of (artificial) boundary condition between the buffer zone and the floodplain. This kind of condition will be discussed in section 3.2.

At the end of K sub-timeloops, we define (h, q) n,m+1 j the 2D water height and discharge per unit width at the lateral boundary vertices j ∈ V b as following:

h n,m+1 j = h n,m,K j , q n,m+1 j i∈L 1D b,j |b ij | = 1 K K-1 k=0 i∈L 1D b,j |b ij |O (1),n,m,k ij , (3.6)
in which O . These values will be used as lateral boundary condition for the black-box (implicit) 2D solver, e.g. Telemac, whose the scheme can be formally written as

V n,m+1 i = V n,m i + ∆t 2D M 2D imp V n,m i , V n,m+1 i ; (h, q) n,m+1 j∈V b , (3.7)
where M 2D imp denotes formally the (implicit) discretization used by the 2D solver. It is well-known that for the shallow-water model, either water height or discharge has to be given if the 2D flow near the bank is sub-critical while both (h, q) n,m+1 j need to be imposed for incoming super-critical flow. In these two cases, Telemac requires also to prescribe the 2D flow direction, namely velocity profile, at the boundary. For this end, one can impose u n,m,K j∈V b , i.e. the 2D flow direction in the buffer zone at the end of these sub-timeloops. Finally for the case of out-leaving super-critical flow, neither water height nor discharge at the boundary is needed.

For the river, finite volume scheme (2.8) can be recast under the form

W n+1 i = W n+1 i - ∆t 1D |I i | b=l,r (O 1D b ) n i , (3.8)
where W n+1 i is the (uncoupled) updated state given by the black-box 1D solver, e.g. Mascaret, (O 1D b ) n i is the rate of overflow at the cell I i and over the lateral boundary b ∈ {l, r}. It is given by

(O 1D b ) n i j∈L 2D b,i |b ij | = 1 M K M -1 m=0 K-1 k=0 j∈L 2D b,i |b ij |O (1,2),n,m,k ij .
(3.9)

A simple transparent boundary condition

To complete the proposed lateral coupling algorithm, a last technical question remains: how wide is the buffer zone? Intuitively, the width of Ω b is related to the traveling distance of water waves, from the lateral boundary into the floodplain Ω 2D , during a 2D timestep ∆t 2D in order to avoid the wave reflecting within the buffer zone. Nevertheless, it is difficult to quantify such distance since, in general, the wave speed depends on time and on the direction of the flow. A more practical strategy proposed here consists in setting an arbitrary width for Ω b , and implementing a reasonable non-reflecting (transparent) boundary condition along the internal (artificial) boundary between Ω b and Ω 2D . This later is nothing but providing the reconstructed states V n,m,k i , i ∈ E ext j in (3.4) when computing the numerical flux for the bank model. To address such an issue, one can use the Riemann invariants associated to the outgoing characteristics as described in the following.

In the present work, the smallest buffer zone is used which consists to define that V b is the set of 2D vertices lying on the lateral boundary b ∈ {l, r}. Let us consider now a face e ji , i ∈ E ext j , i.e. which lies on the artificial boundary between the buffer zone and the floodplain as sketched out in Fig. 4. Regarding numerical flux in (3.4), this comes down to prescribe the reconstructed state V n,m,k j,i , or precisely its projections V L,R on the outward normal direction n ji = (n x , n y ) of e ji , in order to compute the homogeneous numerical flux F b . Indeed, by denoting the projected states

V L,R = R n V n,m,k j,i with R n =   1 0 0 0 n x n y 0 -n y n x   , the 2D numerical flux F b is given by F b (V n,m,k j , V n,m,k i ; n ji ) = R -1 n F 1D b (V L , V R )
with F 1D b the well-known HLLC flux [START_REF] Toro | Shock-Capturing Methods for Free-Surface Shallow Flows[END_REF] in one dimension. Thanks to the simplicity of the hydrostatic reconstruction (2.7), the left-side reconstructed state V n,m,k j can be built by only providing the topography step z jz i , so V L is well computed. Consequently, it remains now to derive the right-side reconstructed state V R . At the first subtimeloop, the available 2D state is set to

V n,m,0 i = V n,m i , i ∈ E ext j , (3.10)
so one can easily deduce V R for k = 0. For k ≥ 1, i.e. when the bank model needs more than one sub-timeloop, and in order to let the waves freely propagate through the face, we seek a right-side state as some small perturbation of the left-side one. As such, we propose to use the Riemann invariants of the shallow-water equations linearized locally around V L := (h, hu, hv) L . In turn, we consider the following system for the first-order perturbation variable ( h, ũ, ṽ):

∂ t   h ũ ṽ  +   u L h L 0 g u L 0 0 0 u L   ∂ x   h ũ ṽ  = 0,
completed with a (partial) initial condition for the associated Riemann problem ( h, ũ, ṽ)(t = 0, x < 0) = 0.

(3.11)

The characteristics dx dt = λ 1,2,3 along which Riemann invariants w 1,2,3 write [23] The reconstruction of V R is formulated as a function of the characteristic wave speed λ 1,2 , as illustrated in Fig. 5 : • Case 1: λ 1 ≤ λ 2 ≤ 0, i.e. all characteristics enter the buffer zone. The right-side state needs to be fully prescribed. To this end, we impose the same state as the one imposed at the first sub-timeloop as in (3.10)

λ 1,2 = u L ∓ gh L , λ 3 = u L , w 1,2 = ũ ∓ h g/h L , w 3 = ṽ. x t 0 λ 1 V L V R λ 2 x t 0 V R λ 1 λ 2 V L x t 0 λ 1 V L V R λ 2
V R = R n V n,m,0 i .
• Case 2: 0 ≤ λ 1 ≤ λ 2 , i.e. all characteristics leave the buffer zone. The right-side state has no incidence on the wave propagation. In order to be consistent with the linearized equations, one can set V R = V L .

• Case 3: λ 1 ≤ 0 ≤ λ 2 , i.e. the second characteristic leaves the buffer zone while the first one enters the buffer zone. From the initial condition (3.11), one has ũ + h g/h L = 0 along dx dt = λ 2 .

(3.12)

Therefore V R needs only to be prescribed partially with water height or velocity. We note also that, on one hand, there might have two faces between C j and C i due to the construction of vertex-centered method cells, on the other hand, the water height remains unchanged when projecting V n,m,k i on these faces. So, let us privilege the water height by defining h R such as

h R = h n,m,0 i , hR = h R -h L , u R = u L + ũR = u L -hR g/h L , (3.13)
which respects (3.12) along the outgoing characteristic. Moreover, it will show in section 3.3 that such definition of h R allows the coupling method to be well-balanced and to preserve the overall positivity of the water height. Finally, the transverse velocity v R is defined depending on the third characteristic

v R = v L if u L ≥ 0, v τ otherwise (3.14)
where v τ stands for the transverse velocity of the projected state R n V n,m,0 i .

Properties of the proposed method

Now let us conclude the present bank model by stating its main properties since they directly result from the construction of the method: Proposition 1. The proposed lateral coupling method using buffer zone by schemes (3.3), (3.7), (3.8) together with exchange source terms (3.5), (3.6) and (3.9) (i) is mass-conservative;

(ii) is well-balanced and preserves the positivity of water height;

(iii) allows large time step with ability to reuse existing 1D and 2D shallow-water solvers.

Proof. The last statement (iii) is evident from setting (3.2) on the timesteps, the generic formalism (3.7)-(3.8) of 1D and 2D models, and thank to the designed transparent boundary condition between the buffer zone and the floodplain.

Let us show (i) after which the coupling is conservative. Indeed, from definition (3.9) the quantity of mass and momentum leaving a 1D control volume I i during one timestep ∆t 1D through the lateral boundary b ∈ {l, r} reads

∆t 1D (O 1D b ) n i |I i | = (M K∆t)(O 1D b ) n i j∈L 2D b,i |b ij | = M -1 m=0 ∆t 2D    1 K K-1 k=0 j∈L 2D b,i |b ij |O (1,2),n,m,k ij    .
Considering the case of sub-critical or incoming supercritical 2D flow and regarding lateral boundary condition (3.6), it could be easy to check that the first component in last equation is nothing but the total discharge that the bank model has injected into the floodplain during M timesteps ∆t 2D . For the case of supercritical out-leaving flow, 2D model needs no prescribed value at the boundary; the resulting overflow quantities are directly computed from the 2D states V n,m j , j ∈ L 2D b,i . But interestingly, this is also what has been done by finite volume scheme (3.3) of the bank model when computing numerical flux and so the exchange source term O n,m,k ij according to (3.5). We notice that for sub-critical flow, exchanged momentum (second component of O n,m,k ij ) might not be exactly conserved due to the fact that only discharge has been imposed in this case. Now we turn to show (ii) which concerns the well-balanced property for lake-at-rest. It is sufficient to show that such steady state is preserved by the bank model. In other words, assume that at time t n,m,k the given data at the bank V n,m,k j , V n,m i∈E j and V n i∈L 1D b,j = V 1D (W n i ) are steadystate solutions (see again Fig. 3), we have to show that V n,m,k+1 j = V n,m,k j . By definition (2.7), the reconstructed states satisfy

u = v = 0, h n i = h n,m,k j = h n,m i = h.
In addition, the reconstructed state V R (Fig. 5) made at artificial boundary is also steady-state solution. Indeed, since the given data corresponds to sub-critical flow, (3.13)-(3.14) lead to

h R = h ⇒ hR = 0 ⇒ u R = v R = 0.
Thanks to the consistency of numerical flux F b , finite volume scheme (3.3) reduces to

V n,m,k+1 j = V n,m,k j - ∆t |C j | i∈E j |e ji | 0 g 2 h 2 n ji + ∆t |C j | i∈L 1D b,j |b ij | 0 g 2 h 2 Γ b = V n,m,k j ,
using the trivial equality

i∈E j |e ji |n ji + i∈L 1D b,j |b ij |(-Γ b ) = 0,
where we recall that, by convention, normal vector Γ b points out into the buffer zone. Finally, one could easily check that the positivity of water height is guaranteed because the water heights involved in finite volume scheme (3.5) are all positive and the fact that numerical flux, e.g. the expected HLLC, is assumed to be positivity preserving. The proof is thus achieved.

Application to the Telemac-Mascaret system

Dynamical coupling with OpenPALM

OpenPALM 2 is an open-source, flexible, and powerful dynamic code coupler that has been developed jointly by CERFACS (Centre Européen de Recherche et Formation Avancée en Calcul Scientifique) and ONERA (Office National d' Études et de Recherches Aérospatiales) since 1998, see [START_REF] Buis | Palm: a computational framework for assembling high-performance computing applications[END_REF]. It was originally designed for Data assimilation (DA) algorithms in operational oceanography forecasting; it has now reached a high degree of maturity and stability, with applications ranging from operational DA (oceanography, atmospheric chemistry, hydrology) to industrially-oriented multi-physics modeling (fluid-structure interactions, combustion-acoustics interactions). OpenPALM allows for the concurrent execution of and the intercommunication between programs coded in various programming languages such as Fortran, C++ and Python. It is defined as dynamic for its ability to deal with situations where the component execution scheduling and the data exchange patterns cannot be entirely defined before execution. OpenPALM provides a straightforward parallel environment based on high-performance implementation of the Message Passing Interface (MPI) standard (i.e. MPICH, OpenMPI, LAM/MPI). This interface can perform both data parallelism (i.e. simultaneous execution on multiple cores of the same code for a unique data set with domain decomposition) and task parallelism (i.e. simultaneous execution on multiples cores of multiple tasks for the same or different datasets). OpenPALM applications are implemented via a graphical user interface called PrePALM where the programmer defines the coupled components as units that are arranged along branches, the communications between these components as well as the general coupling framework including loops, conditional executions, resources management and memory management. Given this layout, OpenPALM builds an executable that is launched with an MPI command on dedicated resources. The PrePALM layout for the lateral coupling algorithm is presented in Fig. 6 and detailed in the following.

Description of the practical implementation

Fig. 6 shows the practical implementation within the OpenPALM framework of the lateral coupling strategy. It is composed of two branches that are executed in parallel when running the coupled simulation (task parallelism). Each branch is dedicated to one of the model: the green one is dedicated to the 1D model Mascaret, while the red one is dedicated to the 2D model Telemac. The routines handling the coupling are integrated within the branch MASCARET. The branches are executed from the top to the base and the units along them (the rectangles on Fig. 6) are executed sequentially. Each unit represents a subroutine and they can themselves be executed in parallel (data parallelism) allowing for running several 2D models in parallel. This is the case of the units of the branch TELEMAC where the OpenPALM variable NB MODELES 2D stands for the number of 2D sub-models used in this study (one for each lateral floodplain). The dashed lines between the units represent OpenPALM communications, i.e. data exchange between units. The green and red shaded areas represent OpenPALM block structures where units within the structure shared specific variables. In particular, at each coupling iteration the 1D and the 2D states are saved within their respective block structure so that it is not necessary to save and read them from a file which ensures the cost-effectiveness of the strategy.

The coupling strategy consists in three steps:

• The initialization consists in reading the user-defined parameters the model parameters, and the 1D and 2D geometries (units " init" in Fig. 6). These variables are allocated and stored in the OpenPALM block structure shared memory;

• The time loop consists in the parallel integration of the 1D and 2D models (units " iter" in Fig. 6). At each coupling step, 1D and 2D data are shared with the coupling module (embedded within the branch MASCARET) which in turns provides the lateral input for the 1D model and the boundary conditions of the 2D sub-models. The data exchanges between the two branches are achieved with OpenPALM communications.

• The finalization step consists in deallocating the variables, closing the output files, saving the 1D and 2D hydraulic states.

The task parallelism as well as using the OpenPALM block structure to manage the variables ensures the numerical efficiency of the implementation of the coupling strategy. The time of execution of the simulation (or wall-clock time) is then approximately that of the largest 2D model.

Validation and first application

Numerical validations aim at highlighting the merits and limitations of the proposed coupling algorithm. This end being devoted to flooding simulation, the validation focuses on the flow in the floodplain and on the lateral discharge between the river channel and the floodplain. The 1D-2D coupled solution is compared to a reference that is either a full 2D solution [START_REF] Goutal | A 2d reconstruction for the transverse coupling of shallow water models[END_REF][START_REF] Morales-Hernández | A conservative strategy to couple 1d and 2d models for shallow water flow simulation[END_REF][START_REF] Morales-Hernández | A riemann coupled edge (rce) 1d-2d finite volume inundation and solute transport model[END_REF][START_REF] Morales-Hernández | Conservative 1d-2d coupled numerical strategies applied to river flooding: The tiber (rome)[END_REF] or an analytical solution when available. Experimental results can also be taken as a reference. In the two latter cases, the validation also takes into account the error related to the numerical modeling choices, on top of the coupling errors.

Here, the coupling algorithm is first validated on a dam-break test case with respect to an analytical solution and to a full 2D solution. It is then validated on a dyke-break test case with respect to experimental results [START_REF] Roger | Experimental and numerical investigations of dike-break induced flows[END_REF] and to a full 2D solution. It is finally validated with respect to the full Telemac2D solution for the historical 1981 flood event over the Garonne River (France).

Dam-break flows in transverse direction

A rectangular channel of 1m width cross-section connects laterally with a rectangular floodplain of 3.5 × 4.0m 2 as illustrated in Fig. 7. The 1D-2D coupling interface is parallel to the channel center line (-2 ≤ x ≤ 2, y = 0). The x and the y axis are set respectively parallel and perpendicular to the hydraulic axis ξ, and centered at the middle of the lateral floodplain. Both channel and floodplain are frictionless and have the same bottom elevation z b = 0m. At the initial condition, the velocity of the flow in the channel is null and the water height is h 0 = 0.5m while the floodplain is empty. Furthermore, open boundary condition are prescribed on the channel. The numerical simulations are carried out over a short period so that the dam-break waves do not reach the left and right boundaries of the channel. Consequently, the prescribed boundary conditions in the flood plain have no influence on the solutions. The 1D and 2D numerical solutions depend on the following numerical settings. A coarse mesh with uniform resolution ∆x 1D = 0.1m is used for the channel, while the floodplain and the full 2D simulation are represented with a finer unstructured triangular mesh of mean resolution ∆x 2D = 0.02m. As such, a 1D cell is adjacent to about five 2D cells at the lateral boundary. The trans-critical kernel [START_REF] Goutal | Mascaret: A 1-d open-source software for flow hydrodynamic and water quality in open channel networks[END_REF] of Mascaret is used to solve the dynamics of the flow in the channel with an explicit finite volume scheme with Roe-type numerical flux to solve (3.8). The timesteps are set to ∆t 1D = ∆t 2D = 0.01s which ensures the CFL condition for the 1D scheme, while preserving the stability of the 2D implicit solver (3.7) of Telemac. It has been observed that the bank model performs about 10 sub-timeloops of the finite volume scheme (3.3) for each timestep; meaning that M = 1 and K 10.

As long as the 2D waves, propagating from the lower corners of the 1D-2D interface, have not yet reached the middle of the floodplain, exact solution for the water height and the velocity along y-axis remains purely 1D and are self-similar (in y/t). It is in agreement with the Stokes solution for the rarefaction wave propagated backward in the channel and downward in the floodplain. For -√ gh 0 ≤ y/t ≤ 2 √ gh 0 , the solution writes:

h(0, y, t) = 1 9g 2 gh 0 - y t 2 , v(0, y, t) = 2 3 gh 0 + y t . (5.1)
In particular, the solution at y = 0 remains constant in time. Therefore, the lateral discharge per width unit (unit-discharge) reaches: Fig. 8 displays the water height along the y-axis in the 1D and 2D domain at t = 0.4s on the left panel and the lateral unit-discharge along time on the right panel. In Fig. 8a, there is a good agreement between the the water heights computed over the 2D domain, from the 1D-2D coupled model, the full 2D model and the Stokes solution (5.1-5.2). Same conclusions are drawn from Fig. 8b for the lateral discharge. The oscillations observed at early stage of the simulation for the full 2D as well as the coupled model are likely to come from the linear approximation of the solutions in Telemac.

hv(0, 0, t) = 8 27 gh 3 0 0.328 m 2 /s. (5.2)
We recall that the 1D model (2.2) considers a constant water level over the cross-section, as one can observe in Fig. 8a. It is why the 1D water level can not be taken as a relevant approximation of the water level at the bank. In other words, using the 1D solution to directly prescribe the boundary condition for the 2D shallow-water equations is not a valid strategy. A proper alternative, proposed here, is to impose the mass flux while maintaining Riemann invariant over the lateral boundary. Finally, it should be noted that the water height near the lateral boundary given by the coupled model remains lower than that of the full 2D solution. This results from the approximation of the transverse velocity close to the bank in the reconstructed 2D state (2.4) which is set to zero, thus leading to the underestimation of the mass flux through the bank parallel to the river.

Dyke-break induced flood wave propagation

This test case is based on a physical model proposed in [START_REF] Roger | Experimental and numerical investigations of dike-break induced flows[END_REF]. It is similar to that presented in 5.1, with a smaller lateral boundary (0.7m) between the 1D and 2D domains, as illustrated in Fig. 9, in order to represent a dike failure. The 1D model boundary conditions are: Q(t) = 0.3m 3 /s at the inlet of the channel and h(t) = 0.5 at the outlet that represents a weir. The 1D initial hydraulic state in the channel is prescribed from the steady state imposed by these boundary conditions (u = 0.6m/s) and the floodplain is supposed to be empty at initial time. The friction coefficient is set to n = 0.015s/m 1/3 for both 1D and 2D domains.

In the dyke-break configuration, the wave front in the floodplain features some deflection from the main y-axis, that was not observed in the dam-break test case. This feature is observed with both the numerical and the physical model. However, it has been showed that the deflection is in general overestimated by the numerical simulation, either full 2D model or higher dimension models including complex turbulent description [START_REF] Roger | Experimental and numerical investigations of dike-break induced flows[END_REF][START_REF] Stilmant | Dike-break induced flows: a simplified model[END_REF][START_REF] Sun | Processes of dike-break induced flows: A combined experimental and numerical model study[END_REF]. This test case is challenging one for 1D-2D coupling validation as it features a complex flow dynamics. The validation is carried out for water height at three locations in the flood plain denoted by P1, P2, P3 in Fig. 9.

Fig. 10 shows the water height and the stream traces at time t = 1s. Full 2D and 1D-2D coupling results are in good agreement and both represent the flow deflection in the floodplain. It should be noted that the deflection is more important in the 1D-2D coupling result. This is consistent with coupling hypothesis made on transverse velocity at the bank. Indeed, the deflection is mainly controlled by transverse velocity along the lateral boundary that is, however, neglected in the 2D reconstruction as written in Eq. (2.4). Furthermore, the full 2D water height varies more significantly in the vicinity of the lateral boundary than that of the 1D-2D coupled solution. This is due to the fact that the local behaviour suffers from 1D approximation that prescribes constant water levels along cross sections of the channel. As already observed and discussed in other works, e.g. [START_REF] Roger | Experimental and numerical investigations of dike-break induced flows[END_REF], the accuracy of the full 2D model is not perfect and the turbulent and/or non-hydrostatic pressure models may need further improvement. The 1D-2D coupling results are closer to the experimental results at P2 and P3, further from the lateral coupling interface than P1 where the water height is significantly under-estimated by the coupling model. Details on the resulting RMSE are reported in Fig. 11b. These conclusions are in agreement with that of the dam-break test case and highlight again the expected limitations due to the transverse velocity reconstruction.

Real applicative case : 1981 Garonne flood event

The 1D-2D lateral coupling is applied to a reach of the Garonne river during the flood event in December 1981 (December 12 th to December 21 st ). This event is ranked among the 10 major flood events since 1875 and was established as a benchmark for hydrodynamic models validation. This case presents with complex bathymetry over 50 km of river between Tonneins (upstream) and La Réole (downstream), as shown in Fig. 12a. A system of longitudinal dykes along the main channel was progressively constructed since the historical flood of 1875 to protect the floodplains, organize submersion and flood retention areas.

Both a 1D Mascaret and a full 2D Telemac models were developed and calibrated over the watershed [START_REF] Besnard | Comparaison de modèles 1d à casiers et 2d pour la modélisation hydraulique d´une plaine d´inondation -cas de la garonne entre tonneins et la réole[END_REF]. In the 1D model, the geometry of the river channel consists in 83 cross-section profiles while the floodplains were simplified by a system of separated storage areas. In this configuration, overflowing from the river to the floodplains and vice-versa are described following with exchange laws. For the 2D model, a triangular unstructured mesh is used, with an increased resolution near the dykes, as shown in Fig. 12b. In the river channel, the mesh is directional, the 2D bathymethy was linearly extrapolated from the 1D cross-section profiles over this mesh. Related studies highlighted that the simulation results depend essentially on a good representation of the longitudinal banks and dykes for both the 1D and the 2D models. We refer to [START_REF] Besnard | Comparaison de modèles 1d à casiers et 2d pour la modélisation hydraulique d´une plaine d´inondation -cas de la garonne entre tonneins et la réole[END_REF] for further details on the model settings and an in-depth discussion.

In the present study with 1D-2D lateral coupling model, the previously 1D and 2D models were modified in order to include the longitudinal dykes in the 1D model geometry and thus reduce the computational cost of the 2D simulation. The geometry of the 1D model was linearly extrapolated from the 83 cross-section profiles in Mascaret to the 2D geometry. The floodplain is divided into 8 sub-domains as shown in Fig. 12c denoted as LFP [1:5] on the left bank of the river and RFP [1:3] on the right bank of the river. A constant resolution ∆x 1D = 100m has been used for the 1D mesh while it has been changed for the 2D model: the mesh size for the main channel does not exceed 80m, that of the dykes is 40m and that of the floodplain is 150m [START_REF] Besnard | Comparaison de modèles 1d à casiers et 2d pour la modélisation hydraulique d´une plaine d´inondation -cas de la garonne entre tonneins et la réole[END_REF]. A common timestep ∆t 1D = ∆t 2D = 10s is used for both 1D (this guarantees the stability of the 1D explicit finite volume solver) and 2D (implicit finite element solver) numerical models.

Numerical simulations over the 10-day flood event in 1981 were performed with the full 2D and the 1D-2D coupled models using common initial and boundary conditions. An observed rating curve is used to describe the downstream boundary conditions at La Réole observing station. An hydrogram is prescribed at the upstream boundary condition at Tonneins; it results from the use of observed water height translated by the observed rating-curve at the Tonneins observing station. Initial conditions for both 1D and 2D equations, in water height and velocities, were established from the non-flooding permanent flow corresponding to the initial discharge Q = 1830m 3 /s of the upstream forcing. The initial water height at the river center line is displayed in Fig. 13a with a blue line for the 1D model and with a red line for the full 2D model. It should thus be noted that the flood plain is dry at the beginning of the simulation. The upstream forcing shown in Fig. 13b reaches main river bed overflow discharge (Q = 2500m 3 /s) at the very beginning of the event, and the bank overflow discharge (Q = 3600m 3 /s) at Day 1. The water flow then occupies the flood plain to reach the flood peak at Day 5 with (Q ≈ 6000m 3 /s) when the recession starts (for approximately 3 days).

The full 2D and the 1D-2D coupling instantaneous water height fields are displayed in Fig. 14a, Fig. 14b and Fig. 14c at Day 2, Day 6 and Day 10 respectively from top to bottom. For each date, the top panel is the full 2D result and the bottom panel is the 1D-2D coupled model result. This illustration shows a good agreement on the floodplain filling dynamics over the flood event, for both filling and emptying stages.

A quantitative assessment is provided in Fig. 15a that displays water height time series for full 2D model and the coupled model at the validation point in LFP3 (see again Fig. 12c). The solutions are in very good agreement during the filling phase. The coupled solution tends to overestimate the flood peak and empties significantly slower than the full 2D. The RMSE over the entire event is 0.197m. A possible reason for the poor 1D-2D coupled solution may be the approximated bathymetry resulting from the linear interpolation of z b at the 83 cross-section profiles, along 50 km of the river banks. This seems to be a very sensible input data of the coupling model. More investigations are needed for forthcoming work.

Finally, it is interesting to study the performance, in term of computational time, of the coupling model for a real case application. It should be noted that the implementation of the coupled solution could be further improved in terms of parallelism; making this solution even more attractive with respect to the full 2D model. For that reason, we here focus on the number of sub-timeloops performed by the bank model during the simulation which can give some idea on the performance of the coupling strategy. Fig. 15b, displays the number of sub-timeloops performed by the bank model over the filling and emptying stages of the flood event. One can find that more than 100 sub-timeloops have been performed during the filling and emptying stages. This can understand that a timestep ∆t ≤ 0.1s has been used when executing the explicit finite volume solver on the buffer zones. Nevertheless, the part of computing time corresponding to the bank model must not be important because the buffer zones have very small size compared to the overall floodplains. Therefore, the proposed approach takes advantage of both the finite volume solver, for estimating the overflow source term, and the implicit solver used on the 2D domain of the floodplains.

Discussion and perspectives

We have proposed in this paper an approach for the lateral coupling of 1D and 2D shallow-water equations dedicated to the modeling of overbank flows between a river channel and its surrounding floodplains. The main novelty of the method consists in setting up for each lateral interface a bank model, i.e. a small part of 2D domain along lateral boundary, namely buffer zone, on which explicit volume solver has been used in order to accurately estimate the exchange source terms of mass and momentum. Performing sub-timeloops only on these (small) buffer zones allows to relax the CFL restriction on the global (large) timesteps imposed by the 1D and 2D models. Furthermore, the present coupling strategy offers the possibility to make use of existing 1D and 2D shallow-water solvers as black-box. The propose approach has been implemented for the open-source industrial hydro-informatic suite Mascaret-Telemac, and the dynamic code-coupler OpenPALM has been used in order to manage data exchanges between the 1D, 2D and the bank models.

The coupling algorithms was first validated over simple configurations such as dam-break and dyke-break flows with respect to analytical solution and experimental data. Numerical results showed that the overflow discharges is well estimated by the coupled model, as expected, while it tends to underestimate the water height near lateral boundaries. Better reconstruction of the transverse velocity might improve the results. Finally, a real application with a 10-days simulation for the 1981 flood event on the Garonne River was presented. Reasonable agreement between the full 2D simulation and the result of the coupled model has been found. In particular, excellent result can be observed for the filling stage, i.e. when the river overflows into floodplain. For the emptying stage, better accurate 1D bathymetry at the river bank might be needed. More investigation and calibration on the coupled model are needed in forthcoming work. The implementation of the coupling strategy has proven to be cost-effective and thanks to the task parallelism provided by the OpenPALM software, the time of execution of the whole simulation is approximately that of the largest of the 2D model. The wall-clock time could further be reduced by running each 2D model on several processors if they are large enough. If the 2D models are smaller than the optimal sub-domain size when running Telemac in parallel (which is roughly 1000 grid points according to Telemac users' feedback, see also [START_REF] Barthélémy | On operational flood forecasting system involving 1d/2d coupled hydraulic model and data assimilation[END_REF]), it is of no use to run them in parallel (which is the case here). Additionally, in the current implementation, all the 2D models are run whether they are wet or not. The total cpu-time could also be reduced by running the 2D models only when they are wet. In this configuration, with no overflowing the computational cost would reduce to that of the 1D model.

Figure 1 :

 1 Figure 1: (a) Plan view for lateral coupling between the river and surrounding floodplains; (ξ, ζ) represents the system of coordinates for the river. (b) Illustration of given data for the Riemann problem at the bank b(ξ).

  1D b ) and (V 2D b , z 2D b ) on each side of the river bank, as illustrated on Fig. 1b. Consider a bank b ∈ {l, r}, one can directly define (V 2D b , z 2D b ) := (V, z b )(t, x, y) with (x, y) ∈ Ω 2D ∩ b, i.e. the 2D hydraulic state and the bathymetry of floodplain closing to the bank. Unfortunately, the state on the side of the river is not yet clearly determined. The 1D hydraulic state is commonly set to z 1D b (ξ) = z b (ξ, b(ξ)), i.e. the river bathymethy at the bank and a reconstruction (ξ, W) → V 1D

Figure 3 :

 3 Figure 3: Schematic of the proposed lateral coupling procedure.

( 1 )

 1 ,n,m,k ij stands for the first component of O n,m,k ij

Figure 4 :

 4 Figure 4: Zoom on the artificial boundary, i.e the boundary of buffer zone Ω b (shadow sub-domain) located inside the floodplain Ω 2D . The red dashed lines stand for the dual mesh on which a vertex-centered finite volume discretization has been used.

Figure 5 :

 5 Figure 5: Riemann problems of the linearized shallow-water equations.

Figure 6 :

 6 Figure 6: Illustration of the implementation of the coupling strategy with the OpenPALM software.

Figure 7 :

 7 Figure 7: Domain setup for dam-break test case: the dam is located along the red dashed line which is also the 1D-2D coupling interface.

Figure 8 :

 8 Figure 8: Results of dam-break test case: analytical solution (black lines), full 2D (red curves) and 1D-2D coupling (blue curves). The initial water height is plotted with a dotted curve.

Figure 9 :

 9 Figure 9: Physical model setup for dike-break test case: the length of dyke failure is 0.7m centered at (x = 0, y = 0). Validation points in the flood plain are located at x = 0 and y(P 1) 0.8, y(P 2) = 1.3, y(P 3) = 2.3.

Figure 10 :

 10 Figure 10: Water heights and stream traces for dyke-break test case at t = 1s given by the full 2D simulation (left plot, channel and floodplain) and the 1D-2D coupling (right plot, only floodplain). The bank model carried out about 10 sub-timeloops for each timestep ∆t 1D = ∆t 2D = 0.01s.

Fig. 11a displays

  Fig.11adisplays the water heights over time for experimental results, full 2D and 1D-2D coupling at validation points P1, P2 and P3 located in the flood plain. As already observed and discussed in other works, e.g.[START_REF] Roger | Experimental and numerical investigations of dike-break induced flows[END_REF], the accuracy of the full 2D model is not perfect and the turbulent and/or non-hydrostatic pressure models may need further improvement. The 1D-2D coupling results are closer to the experimental results at P2 and P3, further from the lateral coupling interface than P1 where the water height is significantly under-estimated by the coupling model. Details on the resulting RMSE are reported in Fig.11b. These conclusions are in agreement with that of the dam-break test case and highlight again the expected limitations due to the transverse velocity reconstruction.

Figure 11 :

 11 Figure 11: (a) Water heights over time at P1, P2 and P3: experiment (black lines with dot symbols), full 2D (red curves) and 1D-2D coupling (blue curves). (b) RMSE of full 2D and 1D-2D coupling models verus experiment

  (a) 50-km reach of the Garonne river from Tonneins (upstream) to La Réole (downstream). (b) Illustration of the full 2D mesh near Sainte Bazeille. Cross-section profiles for 1D model are defined along black lines. (c) 2D domain decomposition into river bed (blue) and surrounding floodplain divided in 8 sub-domains (grey shadowing) denoted by RFP and LPF (Right/Left FloodPlain) used for the lateral coupled model. A validation point is indicated with a red circle near Marmande and Coussin in LFP3 is the location.

Figure 12 :

 12 Figure 12: Domain setup for simulations of the 1981 Garonne flood event.

Figure 13 :

 13 Figure 13: (a) Initial water level along the main channel resulting from permanent flow established with Q = 1830 m 3 /s in both 1D (blue line) and full 2D (red line) models. The 1D bathymetry is plotted with a black line. (b) Observed discharge at Tonneins; main channel overflowing is reached at Q = 2500 m 3 /s, overbank flow and intra-dykes overflow is reached at Q = 3600 m 3 /s.

Figure 14 :

 14 Figure 14: Water height for full 2D (top panels, channel and floodplains) and 1D-2D coupled model (bottom panels, only floodplains).

Figure 15 :

 15 Figure 15: (a) Water height time series for full 2D model and coupled model at the validation point in LFP3. (b) Lateral discharge and corresponding sub-time loops realized by the bank model during the flood event.

http://www.opentelemac.org

http://www.cerfacs.fr/globc/PALM_WEB/

Acknowledgement

This work was initiated during the Cemracs 2016, where the participation of the authors with M. Filonovich was partially supported by the PEPS1 AMIES. In a second time, this work was supported by the H2020 -EoCoE project and the french national program LEFE/INSU. The work of S. Barthélémy was carried out during the tenure of a post-doctoral fellowship supported by the french national department SCHAPI (Service Central d'Hydrométéorologie et d'Aide à la Prévision des Inondations). Sébastien Barthélémy received also support from the Trond Mohn Foundation, under project number BFS2018TMTO1. The authors would like also thank T. Morel for its expertise and help in the use of OpenPALM as well as E. Lepape (SCHAPI) for his kind advises.