Marc Honnorat
email: marc.honnorat@imag.fr>

Joel Marin
email: joel.marin@imag.fr>

Jerome Monnier
email: jerome.monnier@imag.fr>

Xijun Lai
email: laixijin@gmail.com>

Marc Honnorat -Joël

Marin -Jérôme Monnier -Xijun Lai

Joël Marin

Jérôme Monnier

DassFlow v1.0: a variational data assimilation software for 2D river flows

Keywords: Variational data assimilation, parameter identification, river hydraulics, flood Assimilation variationnelle de données, identification de paramètres, hydraulique fluviale, inondations

published or not. The documents may come

Introduction

The St-Venant equations (shallow water) can describe accurately river flows, in the main channel (1D geometry) or in the major bed (flood plain, 2D geometry). Nevertheless, in order to carry out a reliable numerical simulation, models require input parameters such as bed elevation, roughness coefficients in addition to initial and boundary conditions. Generally, these parameters are approximated (eg roughness coefficients), missing (eg inflow boundary conditions or initial conditions) or subjected to uncertainties. Moreover, the governing equations themselves may model ideal configurations, which not necessarily correspond to real configurations. Classicaly, expert-users of computational river flow softwares, see e.g. [?], must perform many simulations (hence spend time) in order to calibrate the numerical model related to a given configuration (trial-error tests). On the other hand, observations from in-situ gauge stations (eg water elevation) or remote sensing observations (satellite images, air photographs, video images) are very important sources of information on the flow state, but they are not qualitatively integrated into the simulation process. To improve the quality of the simulation, data assimilation methods combine in an optimal sense the information from the model and observated data. This allows to identify some parameter values consistent with reality, and /or initial conditions. Variational data assimilation method, see e.g.

[?], [?], is based on the optimal control theory, [?], and it aims to fuse the dynamical model and observations by minimizing a cost function which measures the discrepancy between the simulation results and physical measurements. This method is operational in meterology and since more recently in oceanography. In river hydraulics, variational data assimilation methods have been used for the identification of model parameters in 1D channels, [?, ?], and in 2D, [?, ?, ?, ?, ?, ?, ?, ?]. In real configurations, observations are available only in very small quantities. At best, water levels are measured at very few gauge stations (information very sparse in space); velocity measurements are even more scarce since they require complex human interventions.

DassFlow is a river hydraulics simulation software designed for variational data assimilation on simple hydraulic configurations for research purpose. The model is based on the bidimensional shallow-water equations, solved by the finite volume method using the HLLC approximate Riemann solver. The minimization procedure is based on the L-BFGS algorithm, [?]), which requires the computation of the gradient of the cost function. Since the forward code is written in Fortran 95, the adjoint code is generated by the automatic differentiation tool Tapenade [?] developped by the Tropics project-team at INRIA Sophia-Antipolis. DassFlow software is interfaced with few free and commercial pre and post-processors (SIG tools, mesh generators, visualization tools), which allows to performs computations with real data.

In Section 2, we present the Shallow Water model. In section 3, we present the Variational Data Assimilation method (4D-Var). In Section 4, we detail the finite volume scheme, then the way the adjoint code is automatically generated. In Section 5, we compare some solutions of the forward model with either explicit solutions or real measures (benchmarks INRIA for the forward code). In Section 6, we consider twin experiments which show the capabilities of DassFlow (benchmarks for the data assimilation software). In appendix we detail characteristics boundary conditions and Riemann invariants.

The present research report treats of the mathematical and numerical features of Dassflow only. Concerning the code itself, we refer to the user manual [?] and the developer guide [?].

2 The shallow water model (St Venant)

The equations

The river hydraulics model is based on the bidimensional shallow water equations in their conservative formulation. The state variables are the water depth h and the local discharge q = hu, where u = (u, v) T is the depth-averaged velocity vector. On a domain Ω and for a computational time interval [0, T ], the equations are:

           ∂ t h + div( q ) = 0 in Ω×] 0, T ] ∂ t q + div 1 h q ⊗ q + 1 2 g∇h 2 + gh∇z b + g n 2 q 2 h 7/3 q = 0 in Ω×] 0, T ] h(0) = h 0 , q(0) = q 0 + boundary conditions (see section 2.2) (1)
where g is the magnitude of the gravity, z b the bed elevation, n the Manning roughness coefficient, h 0 and q 0 the initial conditions for the state variables. The quantity c = √ gh denotes the local wave celerity. It is important to notice that if the bed elevation is constant (∇z b = 0) and without the friction term, (1) is an hyperbolic system. Moreover, a transport equation related to a tracer Φ is considered:

∂ t hΦ + div( qΦ ) = s f (2)
with initial and boundary conditions. s f is the tracer source term.

Boundary conditions

We split the boundary Γ of the domain Ω in three parts (see figure 2.1). We denote by:

• Γ in the part where the flow is incoming (type INFLOW)

• Γ out the part where the flow is outgoing (type OUTFLOW)

• Γ wall the part where it is a closed boundary (type WALL)

Figure 2.1: The simulation domain and its boundaries

The user has a large choice of boundary conditions, depending on the boundary condition type.

Wall boundary conditions

On Γ wall , we consider a slip condition:

u • n | Γc (t) = 0 , ∂h ∂n | Γc (t) = 0 ∀ t ∈ ]0, T ]

Inflow conditions

Few boundary conditions are possible at inflow ie on Γ in .

Discharge imposed

We have:

(q • n) | Γin = -q in (t) , ∂h ∂n | Γin (t) = 0 ∀ t ∈ ]0, T ]

INRIA

Incoming characteristics

The characteristic variables of the system are:

w 1 = u • n + g h0 h, w 2 = u • τ and w 3 = u•n-g h0 h, associated to the eigenvalues λ 1 = u 0 •n+c, λ 2 = u 0 •τ and λ 3 = u 0 •n-c respectively. Let us recall that: c = √ gh 0 .
We detail more concerning the characteristics in appendix A. The characteristic w k is incoming if λ k < 0 and outgoing otherwise. The characteristics boundary conditions are implemented as follows:

. For each point of the boundary, we compute the three eigenvalues λ k .

. For each λ k , -if λ k ≥ 0, then w k is incoming; thus it is evaluated from external data, -if λ k < 0, then w k is outgoing; thus it is evaluated from the internal simulation domain.

. Then, the model variables can be retrieved as follows:

-u • n bnd = w1+w3 2 -u × n bnd = w 2 -h bnd = h0 g w1-w3 2

Outflow conditions

Few boundary conditions are possible at outflow ie on Γ out .

Homogeneous Neumann

∂h ∂n | Γout (t) = 0 , ∂q ∂n | Γout (t) = 0 ∀ t ∈ ]0, T ] Water elevation prescribed h | Γout (t) = z out (t) -z b | Γout , ∂(u•n+2c) ∂n | Γout (t) = 0 ∀ t ∈ ]0, T ]
This boundary condition should be applied only in case of sub-critical flow. We refer to appendix B for a justification of the relation

∂(u•n+2c) ∂n | Γout (t) = 0 RR n°0123456789

Incoming characteristics

This case is similar to those related to the inflow boundary Γ in .

Rating curve

A rating curve is a relation between the normal discharge q • n and the water height h : q • n = f (h). Such a relation can be used to specify outflow boundary conditions if combined with the relation:

∂(U+2 √ gh) ∂n = 0.
To this end, we solve the following non-linear system (see figure 2.2):

q • n = f (h) u • n + 2 √ gh = u • n + 2 √ gh int (3)
where u • n + 2 √ gh int is computed from the interior of the domain. 3 Variational data assimilation (4D-var)

The DassFlow software is designed for variational data assimilation (classicaly called 4D-var method). Variational data assimilation is based on the optimal control theory, see [?], [?].

Briefly, this method consists to compute a control vector value k minimizing a cost function which measures the discrepancy between the computed variable and available data.

In our case the control vector k can include the initial condition, boundary conditions, manning coefficient and bed elevation. For example, if inflow discharge is prescribed and water elevation is prescribed, we get: k = (h 0 , q 0 , q in , z out , n, z b ) T

INRIA

If we impose incoming characteristics for both inflow and outflow boundaries, we get:

k = (h 0 , q 0 , [w in 1 , w in 2 , w in 3 ], [w out 1 , w out 2 , w out 3 ], z out , n, z b ) T
We set below the observation function H(k; h, q) and the cost function J(k) = H(k; h, q) which measures the discrepancy between the state of the system (computed variable) and observations. Then the optimal control problem we solve reads:

M in k J(k) (4)
where (h, q) is the solution of the forward model (1).

This optimization problem is solved numerically by a descent algorithm. Thus, we need to compute the gradient of the cost function. This is classicaly done by introducing the adjoint model.

Let us point out that initial conditions, boundary conditions, manning coefficients and bed elevation are only potential control variables. In practice, one manages to identify very few of them only at the same time.

Observations and cost function

Given observations, we define the cost function J which measures the discrepancy between the computed variable (state of the system) and the available data.

In the present version, we assume that water depth and discharge are available at some point and some time instant. In the forthcoming version of DassFlow, we will assume to have some extra observations such as trajectories at surface (lagrangian data).

Given such (eulerian) observations, we define the cost function as the sum of three terms.

• The first term measures the discrepancy between observation data and computed state variable :

J obs (k) = 1 2 T 0 C h h(k, t) -h obs (t) 2 + T 0 C q q(k, t) -q obs (t) 2
where C h and C q are restriction operators.

• The second term measures the discrepancy of the flux: • The third term is a classical regularization term. For example, if we control inflow discharge, we set:

J f lux (k) = 1 2 T 0 G(k, t) -G obs (k, t) 2 RR n°0123456789 with G(k, t) -G obs (k, t) = ∆t |K i | 4 j=1 T -1 ij G 1 (U L , U R ) -G 1 ( ŨL , U R ) with ŨL = [h obs , h obs u, h obs v] T and G 1 (U L , U R )
J p (k) = 1 2 ∂ 2 tt q in 2
Finally, the cost function writes :

J(k) = α obs J obs (k) + α f lux J f lux (k) + α p J p (k) (5)
where α obs , α f lux and α p are scaling coefficients which must be setted by user.

Adjoint model

In order to compute efficiently all partial derivatives of the cost function J(k) with respect to the components of the control vector k, we introduce the adjoint model, [?].

[?].

INRIA

Given the state of the system (h, q), solution of the forward model, we compute (h * , q * ), solution of the adjoint model:

                             ∂ t h * -1 h 2 (q • ∇) q * • q + gh div(q * ) -g q * • ∇z b + 7 3 g n 2 q h 7/3 q • q * = α obs C T h C h h -h obs + α f lux D k (G -G obs ) T (G(k, t) -G obs (k, t)) ∂ t q * + ∇h * + 1 h (q • ∇) q * + 1 h (∇q * ) T q -g n 2 q h 7/3 q * -g n 2 h 7/3 q (q⊗q)q * = α obs C T q C q q -q obs h * (T ) = 0 , q * (T ) = 0 . (6)
Let us point out that the adjoint model is reverse in time, and it includes the observations. Concerning its boundary conditions, for example in the case of inflow discharge prescribed at inflow and water elevation prescribed at outflow, we get:

C.L. : q * | Γin = 0 , (q * • n) | Γ wall = 0 , (q * • τ ) | Γout = 0 , (∂ n h * ) | Γin∪Γ wall = 0 , h * + 2(u • n)(q * • n) | Γout = 0 . (7)
Given the state (h, q) and the adjoint state (h * , q * ), one obtains the following expression of the partial derivatives of the cost function (independent of the derivative of the state):

∂j ∂h 0 (k) = -h * (0) , ∂j ∂q 0 (k) = -q * (0) , ∂j ∂z b (k) = - T 0 div gh(t)q * (t) dt , ∂j ∂n (k) = 2 g n T 0 u(t) h(t) -1 3 u(t), q * (t) Ω dt , ∂j ∂q in (k) = α p ∂ 2 tt T ∂ 2 tt q in -h * | Γq , ∂j ∂z out (k) = (q * • n) c 2 -(u • n) 2 | Γz . (8)
A single integration of the direct model (1) followed by a single integration backward in time of the adjoint model ( 6)-( 7) allow to compute all components of the gradient of the cost function (see figure 3.4).

The optimal control problem (4) is solved using a local descent algorithm, more precisely the L-BFGS algorithm (a quasi-Newton method), implemented in the M1QN3 routine, see [?, ?]. Thus, these partial derivatives are used as input to the minimization algorithm M1QN3. The global optimization process is represented in Fig. 3.4.

In practice, we don't implement the adjoint model presented above but we compute the partial derivatives directly by differentiating the forward code using an automatic differentiation tool, see next section.

Forward code

Adjoint code

Numerical implementation

In this section, we describe on one hand the finite volume scheme solving the forward model and its implementation, on the other hand the automatic differentiation of the forward code in order to obtain the partial derivatives of J. For more details concerning the numerical implementation of the finite volume scheme, one can consult [?], [?].

INRIA 4.1 Forward code and numerical scheme

The forward (or direct) model ( 1) is discretized using by a finite volume method based on the approximate Riemann solver HLLC, see e.g. [?]. The 2D mesh can be a mix of triangular and quadrangular cells.

The finite volume scheme

Equations ( 1) are rewritten as follows :

∂ t U + ∂ x (G(U )) + ∂ y (H(U )) -(S g (U ) + S f (U )) = 0, (x, y) ∈ Ω 2 , t ∈ (0, T ) where U = [h, q x , q y ] T , G(U ) = [q x , q 2 x /h + gh 2 /2, q x q y /h] T H(U ) = [q y , q x q y /h, q 2 y /h + gh 2 /2] T S g (U ) = [0, gh∂ x z b , gH∂ y z b ] T
z b is the topography term, S f is the 2D friction term (Manning law) and g is the magnitude of the gravity.

For the finite volume K i of the mesh, we define the mean value of the state variable U :

U i = 1 |K i | Ki U dΩ , (9)
where |K i | denotes the surface of the cell. We integrate the equation over K i , we use the divergence theorem, we use the rotational invariance property of the equations [?], then we obtain:

Ki ∂ t U dx + Ni j=1 Eij T -1 ij G(T ij U ) ds - Ki (S g (U ) + S f (U )) dx = 0 (10)
where N i denotes the number of faces of the cell K i (3 or 4), E ij is the cell interface and T ij is the rotation of angle θ ij between the normal to E ij and the horizontal plane, Fig. 4.5:

T =   1 0 0 0 cos θ sin θ 0 -sin θ cos θ  
Using the forward Euler scheme for temporal discretization, the equation becomes:

U m+1 i = U m i - ∆t |K i | Ni j=1 T -1 ij G m (T ij U i ) + ∆t S m g (U i ) + ∆t S m+1 f (U i ) (11)
where m = 0, ..., T /∆t is the time index, ∆t is the time step used for the 2D model integration. S m+1 f denotes the friction term treated semi-implicitely in time.

The difficulty is to compute the numerical flux

G m (T ij U i ) = Eij G(T ij U ) ds.
The rotational invariance property of the equations allows to reduce the 2D shallowwater problem to a sum of 1D local Riemann problems. These 1D Riemann problems are solved using the HLLC solver which consists to approximate the 1D flux G, see e.g. [?] and [?]. Moreover, the topography term S g is included in the HLLC solver (see section 4.1.2). Thus, (11) writes :

U m+1 i = U m i - ∆t |K i | Ni j=1 T -1 ij Gm (T ij U i ) + ∆t S m+1 f (U i ) (12)
Only the friction term is treated outside the HLLC numerical flux. It is discretized semi-implicitely in time as follows:

S m+1 (U i ) =     0 -g n 2 u h 4/3 m (q x ) m+1 -g n 2 u h 4/3 m (q y ) m+1    
where n is the manning coefficient.

The local 1D Riemann problem

We detail below the computation of an approximation of the 1D flux G. This is done as follows. 2)

If V = [H, Q n , Q τ ]
T where (τ, n) denotes the tangential and normal unit vectors respectively, we have the local normal 2D flux:

G(V ) = [Q n , Q 2 n /H + gH 2 /2, Q n Q τ H] T
In the homogeneous case (S g = 0), the 2D finite volume solver consists to solve, for each edge, the following local 1D Riemann problem:

∂ t V + ∂ n G(V ) = 0 with V (x, 0) = V L if x n < 0 ; V (x, 0) = V R if x n > 0. The normal flux G(V ) is computed using the HLLC solver, [?].
Finally, we obtain the new solution U m+1 i using (11).

In the non-homogeneous case (S g = 0), the topography term is included in the flux term, using the following modification of G 2 (V ) proposed by LeVeque [?] :

G(V ) = G(V ) +   0 1 2 g(h L + h R )(z bR -z bL ) 0  

Stability condition

Since we use the forward Euler scheme in time, we do not have linear system to solve but the time step must respect a stability condition (CFL type). In the present case, this CFL stability condition is, see [?]:

CFL = ∆t max ( u + c) min(d L,R ) ≤ 1 , (13)
where d L,R is the distance between the cell L and the center of the edge separing it from cell R.

Boundary conditions

In this section, we describe the way we impose the different boundary conditions described in section 2.2.

Let us consider a "boundary cell" L, ie it has no neighboor through its edge i, (see figure 4.6). In the context of a finite volume scheme, we have to evaluate the flux through the edge i. Thus, it becomes natural to define a "ghost-cell" R, for each boundary cell L.

L R

i In order to impose a boundary condition (e.g. inflow dicharge, incoming characteristics etc), we set a particular value to U R = (h R , u R , v R ) T such that the boundary condition is imposed when solving the local Riemann problem through the edge i (ie by computing the numeric flux G(U L , U R )). Let us detail the values U R required for each type of boundary conditions.

Wall conditions On Γ wall , we consider slip type conditions, then:

• Homogeneous Neumann condition on h is imposed by setting h R = h L • Slip condition on u is imposed by setting u R = -u L and v R = v L
Inflow conditions On Γ in , we impose one of the two following boundary conditions.

.

Inflow discharge

The inflow discharge is imposed (unit in m 3 s -1 ). This corresponds to the discharge through the whole boundary and it is distributed in each cell such that the normal velocity u L is constant along the boundary.

-Homogeneous Neumann on h are prescribed by setting h R = h L .

-The inflow discharge is prescribed by setting u R = -q in (t) A n , where A n is the boundary wet section area, and v R = v L .

.

Incoming characteristics

Once incoming and outgoing characteristics w 1 , w 2 and w 3 are computed in the ghost cell as described in appendix A, we set:

-h R = w1-w3 2c -u R = w1+w3 2hR -v R = w 2

INRIA

Outflow conditions On Γ out , we impose one of the four following boundary conditions.

. Water elevation

-We impose h R = z prescribed -(z b ) R , where (z b ) R is the bathymetry in the ghost- cell. -v R = v L -We impose (u R + 2 √ gh R ) = (u L + 2 √ gh L ), that is u R = (u L + 2 √ gh L -2 √ gh R ).
.

Homogeneous Neumann We impose h R = h L , u R = u L and v R = v L .
. Incoming characteristics Once incoming and outgoing characteristics w 1 , w 2 and w 3 are computed in the ghost cell as described in appendix A, we set:

-h R = w1-w3 2c -u R = w1+w3 2hR -v R = w 2
. Rating curve -As described in section 2.2.3, we solve:

q R = f (h R ) u R + 2 √ gh R = u L + 2 √ gh L (14)
-In addition we impose: v R = v L

Adjoint code

In practice, there exists three approaches to compute the adjoint state variable, then the gradient of the cost function.

1. The discretized continuous gradient can be obtained from the adjoint model ( 6)-( 7) discretized using an appropriate numerical scheme which is then implemented.

2. The discrete gradient can be obtained from the adjoint of the direct numerical scheme and its implementation.

3. The "computational gradient" is obtained from the differentiation of the forward code directly.

This last approach ensures a better consistency between the computed cost function (including all types of errors -errors of discretization, rounding errors, iterative algorithms etc-) and its gradient since it is the computed cost function which is differentiated. A large part of this extensive task can be automated using algorithmic differentiation, see [?]. In the case of DassFlow, the direct code is written in Fortran 95 and it is derived using the automatic differentiation tool Tapenade [?].

How to use the adjoint code?

We describe how to define the direct code, then what is the response of the adjoint code automatically generated and finally how to use it.

Let K be the space of control variables and Y the space of the forward code response. In the case of DassFlow, we have :

k = (y 0 , q in , z out ,
dY = ∂M ∂k (k) • dk
The adjoint model is defined as the adjoint operator of the tangent model. This can be represented as follows: ∂M ∂k (k) * : Y ′ -→ K ′ . It takes dY * ∈ Y ′ an input variable and provides the adjoint variable dk * ∈ K ′ at output, see figure 4.9:

dk * = ∂M ∂k (k) * • dY *
Now, let us make the link between the adjoint code and the "computational gradient". • dY * , dk

INRIA dY dk ∂M ∂k (k)
K ′ ×K = dY * , ∂M ∂k • dk Y ′ ×Y (15)
or, using the relations presented above:

dk * , dk K ′ ×K = dY * , dY Y ′ ×Y . (16)
If we set dY * = (0, 1) T and by denoting the perturbation vector dk = (δy 0 , δn, δz b , δq in ) T , we obtain:

0 1 , dy * dJ * Y ′ ×Y =       δy * 0 δn * δz * b δq in * δz * out       ,       δy 0 δn δz b δq in δz out       K ′ ×K
Moreover, we have by definition:

dJ = ∂J ∂y0 (k) • δy 0 + ∂J ∂n (k) • δn + ∂J ∂z b (k) • δz b + ∂J ∂qin (k) • δq in + ∂J ∂zout (k)
• δz out Therefore, the adjoint variable dk * (output of the adjoint code with dY * = (0, 1) T ) corresponds to the partial derivatives of the cost function J:

∂J ∂y0 (k) = y * 0 ∂J ∂n (k) = n * ∂J ∂z b (k) = z * b ∂J ∂qin (k) = q * in ∂J ∂zout (k) = z * out RR n°0123456789
In summary, in order to compute the "computational gradient" (partial derivatives of the cost function J using differentiation of the forward code), first, one runs the direct code then one runs the adjoint code with dY * = (0, 1) T as input.

Automatic differentiation Tapenade [?] is an automatic differentiation tool for Fortran programs. It is devopped by the Tropics team [?] at INRIA Sophia-Antipolis. Tapenade works using source code transformation : it builds the tangent and/or adjoint code automatically from the direct code written in Fortran. One refers also to [?] to get a detailed description of how automatic differentiation works.

Adjoint code validation

We describe below how we check the validity of the adjoint code. Classicaly, we check that it is actually the adjoint of the tangent linear code (scalar product test) and that it computes correctly the partial derivative of the cost function (gradient test).

Scalar product test

The objective of this test is to check if the adjoint code is actually the adjoint of the tangent linear code. In other words, we check the relation ( 15) :

• Given an arbitrary dk ∈ K, we compute using the tangent linear code : dY = ∂M ∂k • dk • Given an arbitrary dY * ∈ Y, we compute using the adjoint code :

dk * = ∂M ∂k *

• dY *

• Then, we compute the following scalar products :

. sp 1 = dY * , dY Y

. sp 2 = dk * , dk K

• And we check if sp 1 = sp 2 or not.

Gradient test

The objective of this test is to check if the adjoint variables dX * computed by the adjoint code correspond to the partial derivatives of the cost function. The Taylor expansion of the cost function j at k for a small perturbation α δk (where α ∈ R + ) writes :

j(k + α δk) = j(k) + α ∂j ∂k (k) • δk + o α δk . ( 17
)

INRIA

We set:

I α = j(k + α δk) -j(k) α ∂j ∂k (k) • δk . ( 18
)
According to ( 17), one must have: lim

α→0 I α = 1.
The gradient test consists to check this property:

• For an arbitrairy k, we compute ∂j ∂k (k) with the adjoint code. • With the direct code, we compute j(k).

• For n = 0, . . . , N :

• We compute α = 2 -n ;
• With the direct code, we compute j(k + α δk) ;

• We compute I α ; • We check if lim α→0 I α = 1 or not.
|I α -1| α #################################### ## TEST DU PRODUIT SCALAIRE ## #################################### Appel du

Water at rest

The first test is the water at rest with a variable topography. We check if the finite volume solver is "well balanced" or not, see e.g. [?]. That means that it must preserves at least water at rest solutions.

Steady solution and Manning-Strickler formula

This test case aims to check the well known Manning-Strickler formula that links water height, discharge and manning coefficient in the case of a steady and uniform state. In the present case (a rectangular and very large channel), the water height h can be computed by: h = q 3/5 n 3/5 L -3/5 i -3/10 (19

)
where q is the water discharge, n the Manning coefficient, L the canal length and i = sin(α).

This test case consists to compare the steady state solution obtained by DassFlow, to the expression (19).

This test case allows to validate or not the balance between the topography term and the friction term.

Dambreak

This test case simulates a dam-break with a flat bottom and without friction. The domain dimension is 50m × 10m, the topography is constant. There is a discontinuity at x = 25m. The initial condition is the following :

h(x, y) = 1 if x ≤ 25 0.1 if x > 25 see fig. 5.12. INRIA and u(x, y) = 2.5 if x ≤ 25 0 if x > 25 x z Figure 5
.12: Initial condition on h for the dambreak test case

The simulation time is 3s and the time step is 0.01. We perform this test case with a regular structured mesh (with n x = 200 and n y = 20) and with an unstructured triangular mesh (based on the same mesh for boundaries). Results are shown in figure 5.13.

Real data test case: Moselle river

In this section, we present a forward run on real data and compare the results with available measures. The study concerns the flood event from Feb.25 1997, 12h00 to Mars.2 1997, 12h00 of Moselle river. The mesh has 2340 cells and 2430 nodes, and it is a mix of triangles and quadrangles, see fig. 5.15. It includes the topography defined at the center of cells and was generated using the ArcGIS software. Time step is dt = 2s. This flood event test case is studied in [?] and data come from [?].

Boundary conditions are imposed using real measures obtained during the flood event. The hydrographs showing the inflow discharge prescribed at inflow boundary and the water elevation prescribed at outflow boundary are shown in figure 5.14.

In order to validate the numerical solution, we compare it to measures available at the EDF gauge station located approximatively in the middle of the main channel, see fig. 6 Benchmarks for the full code (data assimilation)

In this section, we present some twin experiments results. A "twin experiment" means that a first run of the direct model provides observations. Then control variables are changed (we set them as a "first guess") and we run DassFlow in the data assimilation mode in order to try to retreive the control variable values that generated the observations.

Identification of a local topography

Our first twin experiment concerns the identification of the topography in a small scale and academic case. The domain is 30 m long and 4 m large, and the topography is defined by : z b (x, y) = 0.9 exp -1 4 (x -10) 2 exp -(y -1) 2 + 0.7 exp -1 8 (x -20) 2 exp -2(y -3) 2

(20)

The inflow boundary is at x = 0, the outflow boundary at x = 30. Boundaries y = 0 and y = 4 are walls. We use a rectangular structured mesh of dimension 90 × 20. Bed roughness, defined by its Manning coefficient, is uniform (n = 0.025). We impose a constant discharge q in = 8 m 3 /s at x = 0 and a constant water height h out = 1.4 m at x = 30. We obtain a steady state solution after about 80 s of simulation. Figure 6.17 shows the water height of this steady state solution and the topography. From this steady state solution, we extract the forthcoming observations: h obs and u obs every 0.02 s during 20 s on each volume. The objective of this test case is to retreive the topography. The first guess used is a flat bottom. We use the datassim mode, see [?], with the following cost function :

j 1 (z b ) = 1 2 T 0 h(t) -h obs (t) 2 Ω + q(t) -q obs (t) 2 Ω dt , (21)
Figure 6.18 shows the cost function and the norm of its gradient normalized by its initial values, vs iterates (a) and the identified topography (b). We can notice that convergence is obtained and the reference topography is well retreived.

Identification of inflow discharge

We consider a toy test case which includes many features of real cases (eg. Moselle river). The computational domain contains a main channel (river) and floodplains, see figures 6.19 and 6.20).

Again, the present test case is a twin experiment. At the inflow boundary, we set the inflow discharge shown if figure 6.21 (a) simulating a flood event.

Then we perform a forward run to generate observations at points 1 and 2 shown with black stars in figure 6.20(a). Then, we suppose that the inflow discharge is constant (4.95 m 3 s -1 ), and we try to retreive its real value by assimilating observations. We present in Figure 6.21 the identified inflow discharge for different experiments. In Fig. 6.21(a), observations are h and q at each cell and each time step. In Fig. 6.21(b), observations are h at point 1 and (h, q) at point 2, both at each time step. In Fig. 6.21(c), observations are h at point 1 only, but at each time step.

We can notice that the identified inflow discharge is good even with the observation of h at point 1 only. Also, we can notice that the end of the flood event is not well identified. This is the "blind period" phenomena: for example in case (c), the inflow discharge after 270 s can not be identified because the information from the inflow boundary did not reach yet the gauge station.

Identification of Manning coefficients

We consider the same toy test case presented above but we seek to identify manning coefficients in given land-use. Manning coefficients are defined using five different land uses as follows, see figure 6.22:

• In the main channel : n = 0.025 (gravelly main channel)

• Right to the main channel : n = 0.066 (flood plain with bushes)

• Left to the main channel : n = 0.04 (flood plain with little vegetation)

• Near gauging station number 1 : n = 0.03 (pasture, farmland)

• Near outflow boundary condition : n = 0.10 (urban)

Again, this test case is a twin experiment: a forward run generates the observations. For the present test case, the observations are water height at any time at gauge station 1, see figure 6.20(a). As first guess we consider a constant Manning (n = 0.010). The data assimilation process allows to recover perfectly the values used to generate the observations: In a matrix form, this gives:

Area

U t + A 1 U x + A 2 U y = 0 (23) with U = [h, u, v] T , A 1 =   u 0 0 g 0 u 0 0 h 0 0 u 0   , and A 2 =   v 0 0 0 0 v 0 g 0 h 0 v 0   Let n = [n 1 , n 2 ]
T and τ be respectively the normal and the tangent vector to the boundary (see Fig A.23). The matrix If we neglect the variations along x τ , (24) becomes a system of transport equations of w k at speed λ k in the normal direction n. Given an open boundary, w k is incoming if λ k < 0 and outgoing otherwise.

A = n 1 A 1 + n 2 A 2 has 3 eigenvectors: w 1 = u • n + g h0 h, w 2 = u • τ and w 3 = u • n -
Then, the characteristics boundary conditions are implemented as follows:

• For each point of the boundary, we compute the eigenvalues λ 1 , λ 2 and λ 3 .

• For each λ k , -if λ k < 0, then w k is incoming and must be specified from external data.

-if λ k ≥ 0, then w k is outgoing and must be computed from internal data.

INRIA

• Then the state variables can be recovered from characteristics values by:

-u • n bnd = w1+w3 2 -u • τ bnd = w 2 -h bnd = h0 g w1-w3
2

The remaining question is to compute the outgoing characteristics using internal data. Since the finite volume solver is explicit, at current time step n, we need values of the state variable in ghosts cells at time step n -1, see section 4.1.4. This is done using a simple linear extrapolation in the normal direction of the normal, see figure A.23.

B Riemann invariants and boundary conditions

This section details why the relation ∂ ∂n (u • n + 2c) = 0 is imposed if either water elevation or a rating curve is prescribed. We refer to ([?], section 1.3).

Let us consider the 1D shallow water equations with flat topography and without source term: The idea is to transform this equation in the form :

U t + B U x = 0 (
d dt (U + ω) = 0 ( 26
)
where ω is a quantity to be determined. Since gA b is not constant, we have to solve :

dω = gb A dh (27)
In case of a rectangular channel, we obtain: w = 2 √ gh. This means that (U + 2 √ gh) and (U -2 √ gh) are conserved along characteristic curves of direction λ 1 and λ 2 respectively.

For an outflow boundary (u • n > 0), λ 1 > 0 and (U + 2 √ gh) must be computed from the interior of the domain. Since the present finite volume solver is explicit in time, we need values of the state variable on the boundary at time n -1 when computing those of time n. Then, we use a simple extrapolation in the direction of the normal to compute the quantity (U + 2 √ gh) on the boundary. Using a 0-order extrapolation, we find the well-known condition :

∂(U + 2 √ gh) ∂n = 0
Link with the characteristics boundary conditions INRIA

Figure 2

 2 Figure 2.2: Outflow BC: Rating curve and U + 2 √ gh = 0

Figure 3 . 3 :

 33 Figure 3.3: Flux through all edges of cell L

  Figure 3.4: Variational data assimilation process

Figure 4

 4 Figure 4.5: Adjacent cells and use of the rotational invariance property of the equations

Figure 4 . 6 :

 46 Figure 4.6: Internal cell L at boundary and ghost cell R

Figure 4 . 7 :

 47 Figure 4.7: Representation of the direct model. The tangent model becomes ∂M ∂k (k) : K -→ Y. It takes as input variable a perturbation of the control vector dk ∈ K, then it gives the variation dY ∈ Y as output variable, see figure 4.8) :

Figure

  Figure 4.8: Representation of the tangent model.

Figure 4 .

 4 Figure 4.10 (b) shows an example of the scalar product test.

Figure 4 .

 4 Figure 4.10 (a) shows a result of the gradient test. |I α -1| is plotted against α in logarithmic scale. The convergence is good until α > 10 -7 . When α is smaller, the approximation of the partial derivatives is not reliable anymore due to truncation errors, see e.g. [?].

  code linaire tangent... Appel du code adjoint... <Xd

Figure 4 .

 4 Figure 4.10: Adjoint code validation

Figure 5 .

 5 Figure 5.11: Steady state solution for the Manning test case

  Figure 5.13: Dam break test case (all y values are superposed). a) and b) on a regular structured mesh, c) and d) on an unstructured triangular mesh. a) and c) represent h, b) and d) represent u.

  Figure 5.14: Moselle river. a) Inflow discharge imposed. b) Elevation prescribed at outflow boundary.

Figure 5 .

 5 Figure 5.15: Moselle river. Mesh and water depth computed on Feb.26 1997 14h00

Figure 6

 6 Figure6.17: Identification of the topography. Topography and steady state elevation.

  Figure 6.19: Toy test case mesh (a) and bathymetry (b)

  Figure 6.21: In green : reference inflow discharge. In red : identified inflow discharge. (a): Observation of (h, q) everywhere; (b): Observation of h at point 1 and (h, q) at point 2. (c): Observation of h at point 1 only.

g

  h0 h. They are associated to the eigenvaluesλ 1 = u 0 • n + c, λ 2 = u 0 • τ and λ 3 = u 0 • n -c (c = √ gh 0 )respectively. These eigenvectors w 1 , w 2 and w 3 are the so-called characteristic variables.We rewrite (23), using w 1 , w 2 and w 3 :

  Figure A.23: Computation of characteristic variables in a ghost cell using a linear extrapolation of the values in cells (black points).

  25) with U = [h, U ] T , B = U A/b g U , and b(x) is the canal width. The matrix B has 2 eigenvalues λ 1 = U + g A b and λ 2 = U -g A b . The corresponding left eigenvectors are Λ 1 = g, Lagrangian derivative in the direction of the curve of slope λ 1 .

RR n°0123456789

One can notice that the transported quantities obtained are different from those of appendix A. This can be explained by the use of linearized equations in appendix A.

Let us linearize (26) around (h 0 , U 0 ) :

We get:

By setting in (26), we get:

Assuming h ′ and U ′ are small, we find the expression of section A :


Marc Honnorat
email: marc.honnorat@imag.fr

Xijun Lai

Jerome Monnier

Francois-Xavier Le Dimet

J Ér

Ôme Monnier

Franc ¸ois-Xavier Le Dimet

Variational data assimilation for 2D fluvial hydraulics simulations

INTRODUCTION

The numerical simulation of river flows requires a precise modelling of the underlying physics. The bidimensional shallow water equations can describe accurately many free surface hydraulic configurations. However, in order to carry out a realistic simulation of a particular system, the numerical model requires information related to the physical domain. Model inputs such as bed elevation, roughness coefficients, initial and boundary conditions determine the state of the flow and should therefore be determined accurately. Unfortunately, they are usually not well known. Variational data assimilation makes it possible to optimally combine information from the model and observations of the flow in order to identify model parameters that minimize the discrepancy between simulation results and physical measurements. Here, observations of water elevation and velocities are considered.

The method is applied to the identification of boundary conditions and bed elevation in academic test cases and a case with real data. The bidimensional conservative shallow water equations and their discretization using a finite volume scheme are presented in Section 2, the variational data assimilation procedure is described in Section 3. Numerical results on the identification of inlet discharge and bed topography with two academic test cases are presented in Section 3 as well as an application to real data on the Pearl River for the identification of boundary water level. The conclusion is drawn in Section 5.

SHALLOW WATER MODEL

2.1. Shallow Water equations. We use the bidimensional shallow water equations in a conservation formulation. The state variables are the water depth h and the local discharge q = hu, where u is the depth-averaged velocity vector. On a computational domain Ω and for a computational time interval [0, T ], the shallow water equations can be written as

     ∂ t h + div( q ) = 0 in Ω×] 0, T ] ∂ t q + div 1 h q ⊗ q + 1 2 g∇h 2 + gh∇z b + g n 2 q 2 h 7/3 q = 0 in Ω×] 0, T ] h(0) = h 0 , q(0) = q 0 (1)
where g is the magnitude of the gravity, z b the bed elevation, n the Manning coefficient, h 0 and q 0 are the initial conditions for the state variables. In the following, the variable c = √ gh will denote the local wave celerity. To complete this set of equations, one must define boundary conditions. Let Γ be the boundary of the domain Ω, it can be split up in the following way: Γ = Γ q ∪ Γ z ∪ Γ t ∪ Γ w . On boundary Γ q , we prescribe a discharge q and a homogeneous Neumann condition on the water depth. Boundary Γ w corresponds to a wall, we prescribe a slip condition on the velocity and a homogeneous Neumann condition on the water depth. A water elevation zs is prescribed on boundary Γ z and homogeneous Neumann condition for all state variables are prescribed on boundary Γ t :

on Γ q : (q • n) | Γq (t) = -q(t) , ∂h ∂n | Γq (t) = 0 ∀ t ∈ ]0, T ] , (2)
on Γ w :

u • n | Γw (t) = 0 , ∂h ∂n | Γw (t) = 0 ∀ t ∈ ]0, T ] , (3) on
Γ z : h | Γz (t) = zs (t) -z b | Γz , ∂(u•n+2c) ∂n | Γz (t) = 0 ∀ t ∈ ]0, T ] , (4)
on Γ t :

∂h ∂n | Γt (t) = 0 , ∂q ∂n | Γt (t) = 0 ∀ t ∈ ]0, T ] . (5)
The condition (4) on Γ z is valid only for sub-critical flows, i.e. if the local Froude number F r = u c is strictly less than 1. Otherwise, one should prescribe condition (5). 2.2. Finite Volume solver. The bidimensional shallow water equations described above are solved numerically on an unstructured mesh using the finite volume method. The system (1) can be written in a general form as

∂ t U + div F (U) = S(U) , (6)
where U = (h, q) T is the vector of conservative variables, F (U) = G(U), H(U) T the flux vector and S(U) the source term

G(U) =   q x 1 h q 2 x + 1 2 gh 2 1 h q x q y   , H(U) =   q y 1 h q x q y 1 h q 2 y + 1 2 gh 2   , S(U) = 0 -gh∇z b -g n 2 q 2 h 7/3 q .
The computational domain Ω is discretized using triangular or quadrangular cells. We define the mean value of the state variable U on an arbitrary cell K i by

U i = 1 |K i | K i U dΩ , (7)
where |K i | denotes the surface of the cell. By integrating Eq. ( 6) over K i , using the divergence theorem, we obtain

K i ∂ t U dΩ + N i j=1 E ij T -1 ij G(T ij U) ds = K i S(U) dΩ , (8)
where N i denotes the number of faces of the cell K i (3 or 4), E ij is the cell interface (see Fig. 1) and T ij is the 3 × 3 rotation matrix of angle θ ij . The usual flux term derived from the divergence theorem has been replaced by integrals over the cell edges thanks to the rotational invariance property of the shallow water equations (see [Toro, 2001, p. 65]).

Hence, the bidimensional problem actually consists in a sum of one-dimensional Riemann problems that can be solved numerically using a Riemann solver. One can write the following semi-discrete scheme

d dt U i + 1 |K i | N i j=1 T -1 ij G(U L , U R ) = S i , (9)
where G(U L , U R ) is an approximation of the flux through the cell interface E ij . Subscripts L and R denote cells respectively to the left and to the right of the interface.

To compute numerically the discrete flux G(U L , U R ), we use the HLLC approximate Riemann solver [START_REF] Toro | Shock-capturing methods for free-surface shallow flows[END_REF]. This first order scheme handles correctly the transition between sub-critical and super-critical flows, unlike most other first order scheme (see [Zoppou and Roberts, 2003]). The discretization of the bed slope is actually included in the flux term. A forward Euler scheme is used for time discretization. As a result, the following stability condition on the time step ∆t must be satisfied: ∆t

≤ min(d L,R ) max( u +c)
, where d L,R is the distance between the cell center and the center of interface.

VARIATIONAL DATA ASSIMILATION

The state of the flow is solution of the equation system (1) combined with the set of boundary conditions (2-5). Hence, it is determined by the initial conditions (Cauchy problem), by the model parameters (z b and n) and the boundary conditions. The values of these model inputs form the control vector c = (h 0 , q 0 , z b , n, q, zs ) T .

In order to carry out a simulation of a real flow, it is necessary to have a good knowledge of these model inputs. Unfortunately, they are incompletely known in practice, and

R ȳ θij L ≡ K i x y x Eij Figure 1: Two adjacent finite volumes
when an approximation is available, it is often subject to large uncertainties. However, some observations of the flow state may be available, such as water elevations or velocity measurements. These should be in accordance with the simulation results. Therefore, the problem to be addressed consists in identifying the control vector c consistent with both the simulation results and the hydraulic reality represented by observation data. Data assimilation methods make it possible to combine optimally observation data with numerical model in order to improve the simulation results. Variational data assimilation [Le [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]] is based on optimal control theory [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF] and consists in identifying the control vector that minimizes a cost function.

3.1. Cost function. In order to compare the simulation results and the observations of the flow, we introduce a functional

j(c) = T 0 C h h(t) -h obs (t) 2 + C q q(t) -q obs (t) 2 dt , ( 10
)
where c is the control vector, h obs and q obs are observations of the flow, C h and C q are operators that map the state variables to the observations. This cost function measures the discrepancy between the computed state variables and available observations. Its minimization is carried out using a descent algorithm that requires the computation of the gradient of the cost function, i.e. the vector of its partial derivatives with respect to each control variable. We use the quasi-Newton algorithm M1qn3 written by [START_REF] Gilbert | Some numerical experiments with variable storage quasi-Newton algorithms[END_REF]].

3.2. Adjoint model. An efficient computation of the gradient is necessary to carry out an efficient minimization. However, the cost function is not explicitly dependant of the control vector. Indeed, it is an explicit function of the state variable (h, q) which depends implicitly on c through the shallow water equations. Therefore, the gradient of the cost function cannot be computed directly. A direct computation of the gradient using a tangent linear model would require extensive computations. Conversely, the use of an adjoint model makes it possible to compute exactly all components of the gradient with a single backward integration in time of an adjoint model

                         ∂ t h(t) -u • (u • ∇) q + gh div( q) -g q • ∇z b + 7 3 g n 2 u h 4/3 u • q = C T h (C h h(t) -h obs (t)) ∀ t ∈ ]0, T [ ∂ t q(t) + (u • ∇) q + (∇ q)u + ∇ h -g n 2 u h 4/3 q -g n 2 h 4/3 u (u ⊗ u) q = C T q (C q q(t) -q obs (t)) ∀ t ∈ ]0, T [ h(T ) = 0 , q(T ) = 0 , h | Γt = 0 , ∂h ∂n | Γq∪Γw = 0 , q | ∂Ω\Γz = 0 , q • τ | Γz = 0 , h + 2(u • n)( q • n) | Γz = 0 . ( 11
)
The partial derivatives of the cost function are simple functions of the adjoint state variables h and q. For example

∂j ∂h 0 (c) = -h(0) , ∂j ∂q 0 (c) = -q(0) , ∂j ∂z b (c) = - T 0 div gh(t) q(t) dt , ∂j ∂ q (c) = -h | Γq , ∂j ∂ zs (c) = ( q • n) c 2 -(u • n) 2 | Γz .
3.3. Automatic Differentiation. In practice, there are three main methods to obtain an implementation of the adjoint model. The continuous adjoint model ( 11) can be discretized using an appropriate numerical scheme which is then implemented. One can also write the adjoint of the direct numerical scheme and implement it. A better way consists in starting directly from the implementation of the direct model, and writing the adjoint of the code. A large part of this extensive task can be automated using algorithmic differentiation [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF]. Here, we use the automatic differentiation tool Tapenade [Hacoët and Pascual, 2004].

NUMERICAL RESULTS

In order to validate the data assimilation method, we carry out twin experiments. Using a defined set of parameters c ref , some observations are created from the state variable generated by the direct model. They are called (h obs , u obs ). The goal is to identify the value of c ref , from an a priori hypothesis c 0 with data assimilation. 4.1. Validation case: identification of input discharge. Here, we consider the identification of the inlet discharge q, which is a boundary condition in the shallow water equations (1). The domain is a 100 × 8 meters rectangle with an irregular topography defined by: z b (x, y) = 1 2 -x 200 + 1 10 sin xπ 20 cos (y-2)π 6

. The domain is discretized by a regular mesh made up of 100 × 10 cells. The upstream boundary corresponds to the edge x = 0. The Manning coefficient is set to n = 0.025 and a constant discharge q = 5 m 3 /s drives the flow to a steady state, which is used as an initial condition (h ref 0 , q ref 0 ) for the reference flow. A flood is simulated by increasing substantially the discharge at the inlet over a short period of time. From the aforementioned steady state, we carry out a simulation of 80 s with a time step ∆t = 5 100 s and a reference discharge defined by: qref (t) = 5 + 2(t -10) exp -1 100 (t -15) 2 t≥10 s . The corresponding hydrograph is represented with circles in Fig. 2 (a). We can see the propagation of the wave on Fig. 2 (b), where a cut of the free surface in the plane defined by y = 1.6 m is plotted for several simulation steps. The water is flowing in the direction of increasing x. Observations h obs m of the water depth are measured at each time step in a unique observation point located at (x m , y m ). As an initial guess for the discharge, we choose the constant value q0 = 5 m 3 /s. All other inputs are left unchanged. We introduce the following cost function

j 1 (q) = 1 2 T 0 h(x m , y m ; t) -h obs m (t) 2 dt , (12)
where h(x m , y m ; t) denotes the water depth at the point (x m , y m ) at time t. The value of the cost function as well as the norm of its gradient are plotted against the number of iterations of the minimization procedure in Fig. 2 (c). The identified discharge is represented with a continuous line in Fig. 2 (a). We can see that in 35 iterations, the value of the cost function has been divided by 10 8 and that the identified boundary condition perfectely match the reference value.

4.2. Validation case: identification of topography. Variational data assimilation is now applied to the identification of topography, ie. the value of the parameter z b in equations ( 1). We consider a rectangular domain with dimensions 30 × 4 meters and a reference topography defined by

z ref b (x, y) = 0.9 exp -1 4 (x -10) 2 exp -(y -1) 2 + 0.7 exp -1 8 (x -20) 2 exp -2(y -3) 2 .
The domain is discretized by a regular rectangular mesh made of 90 × 20 cells. The inflow boundary is located at x = 0, the outflow at x = 30. The boundaries defined by y = 0 and y = 4 are walls. Thus, the water flows in the direction of increasing x. The Manning coefficient is set to n = 0.025, the water elevation at the outlet is zs = 1.4 m and a constant discharge q = 8 m 3 /s drives the flow to a steady state after about 80s of simulation. Figure 3 (a) shows the reference topography z b and the elevation of the free surface in this configuration. For a simulation period of T = 20 s, observations of water depth observations h obs and velocity u obs are recorded for every cell of the mesh and for every time step.

The aim of the experiment is to retrieve the value of the reference topography z ref b the initial hypothesis of a flat bed z b ≡ 0. To that purpose, the data assimilation scheme is used with the following cost function to be minimized:

j 2 (z b ) = 1 2 T 0 h(t) -h obs (t) 2 Ω + u(t) -u obs (t) 2 Ω dt , (13)
where • Ω denotes the L 2 norm on the whole domain. The value of the cost function and the norm of its gradient, both normalized by their initial values, are plotted against the number of iterations of the minimization process in Fig. 3 (b). One can notice that convergence has been achieved, and that the original shape of the topography has been well retrieved.

4.3. Real case: the Pearl River. In the previous sections, two academic twin experiments demonstrate the good performance of variational data assimilation to identify control parameters using sufficient observations. Now, the method is applied to a real case: the downstream reaches of the Pearl River, Southern China. In the study area (see Fig. 4), the flow is mainly driven by the tidal force. According to the available data, we consider the following data assimilation experiment: given coarse estimates of the initial and boundary conditions, we seek to identify the water level at boundaries BC 1 , BC 2 and BC 6 as well as a better initial state using the observed water level hydrograph h obs at the stations A, B and C. To that purpose, we defined the corresponding cost function to be minimized

j 3 ( zs , h 0 , q 0 ) = 1 2 T 0 3 i=1 |h i (x m , y m ; t) -h obs i (t)| 2 dt (14)
Figure 4: The Pearl river where zs is the water level at the control boundaries BC 1 , BC 2 and BC 6 , h 0 and q 0 denote the initial state variables and h i (t) (i = 1, 2, 3) denotes the water depth at the three observation stations along each cross section at time t. As coarse estimates for water level at the boundaries BC 1 , BC 2 , and BC 6 , we choose a constant value zs = -0.5 m instead of the observation data. Variational data assimilation is then carried out over a period of 36 hours. The decrease of the cost function is rapid in the first three minimization iterations, then slows down with small oscillating gradients (see Fig. 5 (a)). The water elevation at boundary BC 6 is identified successfully, while the identified water level for BC 1 and BC 2 differs significantly from the observations (see Fig. 5 (b)). From the study of the gradient of the cost function, it has been found that its value is one order of magnitude larger at BC 6 than at BC 1 or BC 2 . It means that these two control variables are much less sensitive to the observations at the gauging stations A, B and C than BC 6 . This should be put in parallel with the distance of the measurement locations to the boundaries. Figure 5 (c) shows that when using the identified boundary conditions for a direct simulation, the state variable perfectly match the measurements of water elevations at the three gauging stations.

CONCLUSION

We have presented variational data assimilation for a river model based on a finite volume discretization of the bidimensional shallow water equations. Academic test cases based on twin experiments demonstrate the ability of the method to successfully identify model parameters in various configurations. The inlet discharge in a channel can be identified when an observation of water elevation is available in a single point of the domain. The identification of bed topography is also possible when a dense set of observations of water elevation and velocities is considered. In a configuration with real data on the Pearl River, the method has proved to be able to identify water elevation boundary conditions from sparse measurements of the water level on three cross-sections.

Figure 2 :

 2 Figure 2: Boundary condition identification: configuration

Figure 3 :

 3 Figure 3: Topography identification

-Figure 5 :

 5 Figure 5: Pearl river numerical results

  Synchronous measurements of water level and discharge are partly or fully available with a one hour time interval for the six open boundaries (denoted BC 1 to BC 6 ) as well as for three cross-sections denoted A, B and C. Both discharge and water level observations are available at boundaries BC 4 , BC 5 , BC 6 and cross sections A, B and C, yet only water level has been measured at boundaries BC 1 and BC 2 . The domain is discretized using an hybrid triangular/quadrilateral mesh,


R Hostache

X Lai

J Monnier
email: jerome.monnier@insa-toulouse.fr

C Puech

ASSIMILATION OF ONE SATELLITE SAR IMAGE FOR FLOOD SIMULATIONS. METHOD AND TEST CASE (MOSEL RIVER)

In view of improving numerical flood prediction, a variational data assimilation method (4D-var) applied to a 2D shallow water model and using distributed water level obtained from one Synthetic Aperture Radar (SAR) image is presented. The RADARSAT-1 image leads to water levels with a ±40 cm average vertical uncertainty of a Mosel River flood event (1997, France). Assimilated in the 2D shallow water hydraulic model, these SAR derived spatially distributed water levels prove to be capable of enhancing model calibration. Indeed, the assimilation process can identify some optimal Manning friction coefficients. Moreover, used as a guide for sensitivity analysis, remote sensing water levels allow also in identifying some areas in the floodplain and the channel where Manning friction coefficients are homogeneous. This allows basing the spatial segmentation of roughness coefficient on floodplain hydraulic functioning.

Introduction

High spatial resolution SAR spaceborne images allow the estimation of distributed water levels in floodplains with reasonable uncertainty by merging SAR derived flood extent limits with a high-resolution highaccuracy Digital Elevation Model (DEM), see Ref. 1. Furthermore, to be reliable, hydraulic models have to be constrained by using various observed data sets. The model calibration consists in forcing the model to provide outputs as close as possible to observed data by searching optimal values of its parameters. A "hand" calibration of Manning coefficients is often done through trial tests with the use of point observations, such as recorded hydrographs at stream gauges. The variational data assimilation 49 method based on the optimal control theory of partial differential equations (also called 4D-var method) offers a special powerful tool to fuse in an optimal sense measurements (observations) and the mathematical model. In river hydraulics, variational data assimilation methods have been used successfully for shallow water models. 2-10 The present study aims at investigating whether the variational assimilation of water levels derived from a flood SAR image could help to enhance the calibration of flood inundation models. This would help getting benefits of both variational assimilation and recently developed remote sensing methods, improving flood model calibration.

Flood Description and Data Extraction

The area of interest includes a 28 km reach of the Mosel River between Uckange (France) and Perl (Germany). In this area, the Mosel River meanders in a flat plain having an average width of 3 km and a mean slope of 0.05%. It is worth noting the presence of a narrow valley at downflow. The latter behave as a bottleneck during flood events causing upstream water retention area. The propagation velocity of the flood peak in the study area is low, around 2 km h -1 . The peak discharge recorded at Uckange city stream gauge (upstream boundary) was around 1,450 m 3 s -1 , corresponding to a 4-5 year time return period.

As hydrometric data, discharge hydrographs were available at three stream gauges located in the study area. These hydrographs are shown in Fig. 1. Recorded discharge hydrographs araise from calculation using observed water stage hydrographs and rating curves: relationship between discharge and water level computed using discharge in situ measurements. Considering that higher magnitude discharges have been in situ measured in Uckange than in Perl for the rating curve computation, the Uckange hydrograph has been assumed more reliable. As a consequence, only the discharge hydrographs in Uckange and EDF stream gauges are used as ground truth information in this study. Consequently, available ground observations are fairly limited. The time series data of discharges at the EDF gauge station are only available at the beginning and the end of the flood period; the measurements are lacking during the high flood stage because of sensor disability (see Fig. 2).

The SAR image (Systematic Map Image product), amplitude coded, has a pixel spacing of 12.5 m, resulting from the sampling of a complex Furthermore, five air photographs at 1:15,000 scale, flown by the French National Institute of Geography (IGN) in 1999 out of flooding period, have been acquired. They have been used to digitize land use information: buildings, urban areas and sparse habitat, and high vegetation, forests, sparse trees, hedges, etc. Indeed trees and buildings perturb the backscattering signal, so we choiced to mask these areas before processing: the whole process runs using only reliable areas.

The topographic and bathymetric raw data have been provided respectively by the North-Eastern French Navigation Services (SNNE) and the DIREN Lorraine as 3D points and 3D lines -calculated by photogrammetry using air photographs at 1:8,000 scale and SONAR sounding. The average altimetric uncertainty on the raw data is about 25 cm (DIREN information) for the floodplain and 1 cm (SNNE information) for the channel. Using a linear interpolation between points and lines, a Triangular Irregular Network (TIN) Digital Elevation Model (DEM) has been generated and then converted to a RASTER DEM (pixel size: 7 m), easily superimposable with image data. Hereafter, DEM acronym will be assigned to the topographic RASTER data. Moreover steep banks, potentially badly represented, are eliminated from the final process.

Water Level Estimation Using SAR Image

Based on the method developed in Ref. 11 providing water level estimates with a ±18 cm mean uncertainty using flood aerial photographs, the water level estimation method used in this study is more detailed in Ref. 12 (see also Ref. 13).

Using the RADARSAT image, this approach provides spatially distributed water level estimates within a ±40 cm mean uncertainty. It can seem strange to reach such accuracy with large pixels (25 m). Nevertheless, more than to pixel size, the accuracy in this process relates to DEM accuracy and the fact that the vertical estimates concern only places with very low slopes and no perturbing items (trees, houses, etc.).

Flood extent extraction

To deal with the radiometric uncertainty, the discrimination of flooded and non-flooded pixels is done using two threshold values. The first threshold value, T min , aims at detecting only pixels that correspond to water bodies. T min has been determined as the minimum radiometric value of nonflooded pixels inside grassland areas (outside the floodplain and outside the permanent water surfaces). 12 The second threshold value, T max , aims at detecting all flooded areas, at the risk of detecting in addition nonflooded areas that have a similar radiometric value to the flooded one. T max has been determined as the maximum radiometric value of water bodies outside the flooded area, using the SAR image pixels located inside the Mirgenbach lake. The thresholding of the SAR image using T min and T max provide a flood extent map with fuzzy limits, coded as follows (see Ref. 12), depending on the intensity I of the SAR image pixels: 0, non-flooded (I > T max ); 1, flooded (I < T min ), 2, fuzzy limit (≈ potentially flooded) (T min ≤ I ≤ T max ).

The innovative point of the SAR image processing is the analysis of the relevance of the remote sensing-derived flood extent limits for hydraulic purpose, and especially for water level estimation. To estimate water levels, the flood extent limits are merged with the underlying DEM. Using such a merging, any erroneous flood extent limit will lead to errors in water level estimation. Consequently, flood extent limits prone to error have to be identified before the merging. Errors in the flood extent limits are mainly due to emerging objects such as building and high vegetation 12 that may mask water. To treat this potential errors, it has been chosen to remove all SAR derived flood extent limits located in habitat or vegetation areas. Considering that radiometric and spatial uncertainties have been taken into account, the "enlarged" relevant limits are then assumed to include the real flood extent limits (Hyp. 1 ). These are shaped as small patches which are sparsely distributed along the floodplain (Fig. 3).

Preliminary water level estimation

The second part of the process estimates one range of possible water levels IWLE = [WL min ; WL max ] for each relevant patch. To do so, the maximum and the minimum elevation values are first extracted inside each relevant patch using the DEM Z values. Next, the DEM altimetric uncertainty (uncert DEM ) is taken into account by being, respectively, added/subtracted to the maximum/minimum values extracted previously for each relevant patch: Since the DEM altimetric uncertainty is taken into account, Hyp. 1 allows to assume that each range of water level estimation -IWLE -includes the real water level (Hyp. 2 ).

IWLE = [min(Z patch ) -uncert DEM ; max(Z patch ) + uncert DEM ].

Final water level estimation

The last part of the process uses hydraulic rules to constrain the water level estimates and reduce their uncertainty, as proposed firstly by Refs. 11 and 14. In a floodplain, hydraulic laws manage the flow, so that water levels must follow an hydraulic logic: hydraulic energy decreases from upstream to downstream. With low flow velocity, like in the Mosel floodplain, this hydraulic rule can be simplified into a decrease of water level in the flow direction (Hyp. 3 ). To apply Hyp. 3 on the IW LE intervals, flow directions between patches (locations of the water levels) have to be determined. As proposed in Ref. 12, some flow directions between patches have been determined using the shape of the SAR derived flood extent and the lines perpendicular to the elevation contour lines, oriented from the highest to the lowest elevation, called steeper lines hereafter. Furthermore, the main flow directions are assumed to be convergent toward the river channel, following the steeper lines. The steeper lines around relevant patches have been determined using the contour lines derived from the DEM. Using this criterion, some up-/down-stream relationship between patches have been determined. At the floodplain scale, these relationships constitute an hydraulic hierarchy of the relevant patches. Consequently, according to Hyp. 3, the water level must decrease from the patch A to the patch B if A is upstream of B. Due to Hyp. 2 this induce the following constraints:

N max (B) ≤ N max (A) (constraint on the maxima), N min (A) ≥ N min (B) (constraint on the minima).
To apply these constraints, the algorithm that has been developed, is flow oriented and impose a decrease on the maxima from upstream to downstream, and vice et versa, an increase on the minima from downstream to upstream. This finally provides intervals of constrained water level estimation IWLE = [WL min ; WL max ], with a half mean range of about ±40 cm.

As a consequence, the method allows the definition of a set of distributed water levels across the floodplain at the satellite overpass time. These remote sensing-derived water levels are then used as new observation for the hydraulic modeling.

As a matter of fact, remote sensing derived data are assumed to be an observation of water depth along the floodplain. The single image is extracted from the satellite image at 6:00 am, Feb 28, 1997. The locations where the water level values have been extracted are shown in Fig. 3.

Mathematical Model with Variational Data

Assimilation (4D-var)

For more details, we refer to Refs. 13 and 15. Numerical computations are performed by using the DassFlow software. The full inverse model (with variational data assimilation) includes the forward model (2D SWEs), the adjoint model and the minimization algorithm. The 2D-SWEs is solved using a finite volume method on unstructured meshes. The adjoint model is directly derived from the source codes of forward model by automatic differentiation. The cost function J that has to be minimized is defined in detail below. As a matter of fact, the optimization problem to be solved is: min p J(p), where p is the control variable. Here, the latter is the Manning coefficient (spatially distributed). We refer to Refs. 13 and 15 for more details.

Preliminary forward run

As a preliminary study, the Manning coefficients are set to an a priori "reasonable" constant value (equal to 0.025). Then, the simulation results are compared with the available measurements, namely the water stage hydrograph at the middle gauge station (EDF) and the SAR image derived water levels in the floodplain at the satellite passover time.

The computed water level at EDF gauge station is very close to measurements; while in floodplain at image time, the computed water level fits much less with available data. Nevertheless, this preliminary forward run shows that the numerical model can reproduce the flow; and the current solution will be the first guess solution (see next section).

Let us point out that the Manning value imposed (0.025) is the best value we found by hand, i.e., using a classical trial-error approach. Furthermore, we did not improve significantly the result if using the landuse information or not (and imposing classical Manning values related to the land-use considered). We refer to Ref. 13 for more details.

Cost function

The cost function J to be minimized (see Section 4.1) contains three terms: J obs (p)=J obs + αJ flux + βJ reg (α and β are weight coefficients to be set).

(1) J obs corresponds to the discrepancy between the observations and the computed flow state. (2) J flux corresponds to the discrepancy of net mass flux.

(3) J reg is a regularization term for smoothing time-dependent control variables p.

The term (J obs + αJ flux ) is defined as follows:

J obs + αJ flux = 1 2σ 2 z m z Timag (x m ,y m ) -z obs Timag (x m ,y m ) 2 1+α (u, v)(x m ,y m ) 2 + γ 1 2σ 2 Q T1 0 (Q(t) -Q obs (t)) 2 dt + T T -T2 (Q(t) -Q obs (t)) 2 dt
where σ z and σ Q are the standard deviations of observations, T imag is the image time, (x m ,y m ) is the position of water level measurements, γ 1 is a weighting coefficient, T 1 and T 2 are time periods of the discharge measurement available at EDF station (see Fig. 2), T is the assimilation time period.

Calibration of Spatially Distributed Manning Friction Coefficients (Land-Cover Based Spatial Distribution)

The Manning friction coefficients n, that represents the resistance to the flow in channels and floodplains, are empirical. In fact we must point out that (n) does not really exists, since it is strongly scale dependent, because it integrates all the friction processes at all scale. This section aims at showing the capability of SAR derived spatially distributed water levels to enhance the Manning coefficient calibration, in comparison to a "hand" calibration using trial-error tests. Thus, the data assimilation method presented here uses "orthogonal data": water heights derived from water limits versus classical water depths in gauge station, see Fig. 2. In this section, the spatial distribution of Manning friction coefficient is based on land-cover classes as classically done. Various cases have been considered depending on the total number of classes taken into account (between 1 and 10). The domain is decomposed by using land-cover classes, one Manning coefficient value n being set by land-cover class. Four various decompositions have been investigated, depending on the total class number:

(a) A lumped n value: one constant value all over the domain. (b) Three classes (distributed n) consisting of the main channel and the left and right floodplains. In all numerical experiments, the assimilation period is 66 h from 12:00, Feb 25th (flood event starting time) to 6:00 am, Feb 28, 1997 (SAR image acquisition time).

Identification of Manning coefficients using synthetic data (twin experiments). We performed identification experiments based on synthetic data, see Ref. 13 for more details. They showed that the SAR derived spatially distributed water levels and the stream gauge in situ measurements can identify properly the spatial distributed Manning coefficients, at least if three land-cover classes only are considered and if the mathematical model is perfect (i.e. observations corresponds to the modeled flow) and no error is introduced in measurements. These last two features are very important. Next, we use the real data and we will notice that the (real) identification problem is much more difficult to solve numerically, even with three landcover classes only.

Identification of Manning coefficients using real data. In all cases (depending on the land-cover number), the computed value of the main channel is the same whatever the initial guess n 0 (e.g. 0.033, 0.020, 0.025 or 0.050) and the total class number. On the contrary, n in other land-cover classes are almost not modified during the optimization process. This result is consistent with the forthcoming sensitivity analysis (see next section), runs which shows that the total cost function is almost insensitive to the n values outside the main channel.

Using real data, sensitivities can become different, and the identification problems become much more difficult. In the present case, the few sensitivity of the model to the Manning coefficient outside the channel is mainly due to the rather low time return period of the investigated flood event (4-5 years). The results show that the water levels in the floodplain is more influenced by the water level inside the channel (and thus by the channel Manning coefficient) instead of the floodplain Manning coefficient.

In Figs. 5 and6, we plot the difference between the computed and the observed water levels where image information is available and when considering 1, 3 and 10 land-cover classes. Figure 7 shows the comparison between simulated and observed discharge hydrographs at the EDF middle gauge station (in function of time). These results show that after the identification of an "optimal" set of Manning coefficient values, the computed flow state is much closer to the observations than that computed using the forward hand-calibrated model (i.e. with n values resulting of the trial-error tests-"hand" calibration).

Sensitivity Analysis and Manning Decomposition

The previous section shows that our variational data assimilation process is capable of enhancing significantly the calibration of the distributed Manning friction coefficients, and consequently the accuracy of the numerical model. Nevertheless, all computations have been done with an a priori: the Manning coefficients have been defined spatially distributed according to given land-use classes (between 1 and 10). This a priori is rather traditional. Nevertheless, does this necessary lead to the most accurate numerical model? In the current section, we aim at answering to this question. To meet this aim, a sensitivity analysis is performed without any a priori for the Manning spatial distribution (i.e. any land-use is defined).

One run of the forward model plus the adjoint model, without running the full optimization process, gives the gradient of cost function J with respect to the Manning coefficients. These gradient values represent the local sensitivity of the cost function with respect to the Manning coefficients. They help one to understand which Manning area is the most important to calibrate.

First, we perform sensitivity analysis in the case the Manning is decomposed by the 10 land-uses defined previously, see Fig. 4. We show the 10 spatially distributed gradient values in Table 1. This sensitivity analysis suggests that the most important Manning value to focus on is in the main channel, than much less important are the vegetation area, the bridge (in main channel), gravel area and grassland area. Others land-use values are negligible.

This sensitivity analysis result is consistent with the calibration process presented in previous section: the optimization algorithm calibrate essentially the main channel value. No a priori land-use case

Next, a sensitivity analysis with no a priori on the Manning decomposition has been performed: we define one Manning value for each mesh cell, and we compute the gradient value for each cell. In particular in the main channel, there are as many potential values as cells. The computed spatially distributed gradient values are shown in Fig. 8. The present sensitivity analysis suggests defining few Manning areas inside the main channel (at least four in the present case). More generally, this suggests to define the Manning areas not upon the land-use (or land-cover) but upon "sensitive areas". Furthermore, the combination of the satellite image and runs of sensitivity analysis indicates the most important Manning areas we need to calibrate (including inside the main channel).

Conclusion

By using a variational data assimilation approach (4D-var), we have investigated potential contributions of SAR derived spatially distributed water levels for the identification of time-independent parameters (Manning coefficients) in a shallow-water flood model.

The spatially distributed water levels have been derived from a SAR image by employing the method developed in Ref. 12. They have been obtained with a ±40 cm mean uncertainty, using a RADARSAT-1 image of the 1997 flood event of the Mosel river. This has been possible by using both an analysis of the relevance of SAR derived flood extent limits for hydraulic purposes, and a merge between the relevant limits and a highly resolution DEM, under hydraulic coherence constraints inspired from Refs. 11 and 14. Such a water level estimation provides spatially distributed information at the time of a satellite overpass while classical in situ measurements are punctually available.

Numerical experiments conducted in this study indicate that a rather dense information in space is of great benefit for the identification of unknown parameters (Manning friction coefficient). Indeed, the assimilation of the SAR derived water levels, in addition to an incomplete discharge hydrograph, proves to be capable of identifying Manning friction coefficients, while the ground data alone does not allow such an identification. Furthermore, a sensitivity analysis conducted by using the SAR derived water levels, shows that a spatial distribution of the friction coefficient based on land-cover may not necessarily lead to the better model results. Indeed, these water levels, used as a guide in the sensitivity analysis, can define areas of Manning friction homogeneity, without apparent link with land-cover. Such sensitivity analysis may finally base the Manning parameter spatial distribution more on the model hydraulic functioning, than on the land-cover.

In a near future, with the launch of new radar satellites with better spatial and radiometric resolutions and more suitable wavelength, the uncertainties of water levels estimates will presumably be further reduced.

Fig. 1 .

 1 Fig. 1. Recorded hydrographs during the flood event at three different locations.

Fig. 2 .

 2 Fig. 2. Observations available: discharge at low water level at EDF gauge station and one (partial) image.

Fig. 3 .

 3 Fig. 3. Spatial distribution of water level values available after image post-treatment.

6

 6

4. 1 .

 1 Mathematical modelThe flood flow is modeled by the two-dimensional shallow water equations (2D-SWEs), see Refs. 13 and 15 for more details. The unknowns are h the water depth, Q x = hu and Q y = hv the unit discharge in the x-and y-directions ((u, v) is the velocity). The friction slopes are evaluated using the Manning formula.

Fig. 4 .

 4 Fig. 4. Land cover: 10 classes (Case d) in the text.

Fig. 5 .

 5 Fig. 5. Water levels (m) at image blocks, i.e. where image information is available (see Fig. 3). Vertical bars correspond to the measures with estimated uncertainties. The four curves correspond to the preliminary forward run and the three calibrated model responses (depending on the land classes number considered; VDA-nXX = 1, 3 and 10 land classes) (Recall: VDA = Variational Data Assimilation = 4D-VAR algorithm; direct = forward model without VDA ("hand" calibrated model).

Fig. 6 .

 6 Fig. 6. Difference of water depth (m) at image instant (6:00, Feb 28, 1997) between observations and computed values using the calibrated model (Manning is calibrated in three land classes, Case b).

Fig. 7 .

 7 Fig. 7. Measured discharge hydrographs (m 3 s -1 ) at the middle gauge station EDF, and computed one using the calibrated model (Manning is calibrated in three land classes, Case b; at downstream, elevation is imposed).

Fig. 8 .

 8 Fig. 8. Sensitivity analysis. Gradient value: cost function with respect to the Manning coefficient in each finite volume cell (there is no a priori land-use).

Table 1 .

 1 Sensitivity analysis. Gradient values for 10 land-uses (at Manning value constant everywhere n =0 .033).

	Channel	Vegetation	Bridge	Gravel	Grassland
	881.54 e + 6	43.75 e + 6	17.76 e + 6	3.64 e + 6	-1.06 e + 6
	Left	Right	Snaked land	Urban	Downstream channel
	0.07 e + 6	0	0	0	0

R. Hostache et al.

The capability of SAR derived water levels may be enhanced to help the identification of model parameters using variational assimilation.