Benoit Combès

Dominique Heitz
email: dominique.heitz@irstea.fr

Anthony Guibert

Etienne Mémin
email: etienne.memin@inria.fr

A particle filter to reconstruct a free-surface flow from a depth camera

Keywords: Data-assimilation, Free-surface flow, Particle filter, Ensemble Kalman filter, Depth sensor, Shallow-water equations

We investigate the combined use of a Kinect depth sensor and of a stochastic data assimilation method to recover free-surface flows. More specifically, we use a Weighted ensemble Kalman filter method to reconstruct the complete state of free-surface flows from a sequence of depth images only. This particle filter accounts for model and observations errors. This data assimilation scheme is enhanced with the use of two observations instead of one classically. We evaluate the developed approach on two numerical test cases: a collapse of a water column as a toy-example and a flow in an suddenly expanding flume as a more realistic flow. The robustness of the method to depth data errors and also to initial and inflow conditions is considered. We illustrate the interest of using two observations instead of one observation into the correction step, especially for unknown inflow boundary conditions. Then, the performance of the Kinect sensor to capture temporal sequences of depth observations is investigated. Finally, the efficiency of the algorithm is qualified for a wave in a real rectangular flat bottom tank. It is shown that for basic initial conditions, the particle filter rapidly and remarkably reconstructs velocity and height of the free surface flow based on noisy measurements of the elevation alone.

Introduction

The circulation of water in open channel (or free-surface flow) is a key variable in many problems involved in hydrology. Most often, it is not possible to observe the complete free-surface flow (i.e. time and space dependent velocity and elevation) but only temporal sequences of sparse measurements coming from e.g. drifters, altimeters or satellites. Although insufficient in complex situations, these observations can be used as the starting ingredients to obtain a more complete characterization of the flow under study. A common way to proceed consists in coupling the observations with the conservation laws of the flow using a data-assimilation (DA) method. In geophysics, DA has a long history of developments and applications. In fluid dynamics, DA was first used to provide low order representations of the flow (D 'Adamo et al 2007[START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF]. Then, for the first time, using variational approach and particle filter, turbulent flow image sequences were coupled with direct numerical simulations [START_REF] Papadakis | Variational assimilation of fluid motion from image sequences[END_REF] and with vortex particle simulations [START_REF] Cuzol | A stochastic filtering technique for fluid flows velocity fields tracking[END_REF], respectively. The reader can refer to [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF] for a synopsis of these attempts. More recently, DA has been carried out either with sequential approaches [START_REF] Colburn | State estimation in wall-bounded flow systems. Part 3. The ensemble Kalman filter[END_REF][START_REF] Suzuki | Reduced-order Kalman-filtered hybrid simulation combining particle tracking velocimetry and direct numerical simulation[END_REF] or with variational approaches [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF][START_REF] Foures | A data-assimilation method for Reynolds-averaged NavierStokes-driven mean flow reconstruction[END_REF][START_REF] Mons | Is isotropic turbulence decay governed by asymptotic behavior of large scales? An eddy-damped quasi-normal Markovian-based data assimilation study[END_REF] or even through the interesting combination of both schemes [START_REF] Yang | Enhanced ensemble-based 4DVar scheme for image assimilation[END_REF].

Data-assimilation using noise-free dynamic model

In the context of free-surface flows, many methods propose to tackle the dataassimilation problem by considering that we know the dynamic model describing the evolution of the flow under study up to some parameters such as the bed roughness, the inflow velocity, etc. The goal of the DA method is then to estimate these unknown parameters of the model, fitting as best as possible the observations. In addition, to provide an estimate of the free-surface flow, such an approach is able to handle more specific applications such as: the estimation of equivalent topography (bed geometry and roughness) to fit the parameters of a given dynamic model to a particular river then allowing to perform suitable numerical simulations [START_REF] Roux | Use of parameter optimization to estimate a flood wave: Potential applications to remote sensing of rivers[END_REF]; the forecasting of a given situation by estimating the current flow and the parameters of the model. In the literature, it has been many times shown that such methods successfully perform on synthetic and sometimes real specific application cases (e.g. [START_REF] Hostache | Assimilation of spatially distributed water levels into a shallow-water flood model. Part II: Use of a remote sensing image of Mosel River[END_REF][START_REF] Honnorat | Identification of equivalent topography in an open channel flow using Lagrangian data assimilation Computing and Visualization[END_REF]. As for now, most of these works focused on the estimation of the initial condition of the flow [START_REF] Tinka | Quadratic Programming based data assimilation with passive drifting sensors for shallow water flows in[END_REF][START_REF] Titaud | Assimilation of Image Sequences in Numerical Models[END_REF][START_REF] Lai | Assimilation of spatially distributed water levels into a shallow-water flood model. Part I: Mathematical method and test case[END_REF], the time dependent flow at the open boundaries [START_REF] Strub | Inverse estimation of open boundary conditions in tidal channels[END_REF][START_REF] Bélanger | Data assimilation (4D-VAR) to forecast flood in shallow-waters with sediment erosion[END_REF][START_REF] Castaings | Automatic differentiation: a tool for variational data assimilation and adjoint sensitivity analysis for flood modeling[END_REF][START_REF] Lai | Assimilation of spatially distributed water levels into a shallow-water flood model. Part I: Mathematical method and test case[END_REF][START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations Computing and Visualization in[END_REF], the bed roughness coefficients [START_REF] Ding | Identification of Manning's Roughness Coefficients in Shallow Water Flows[END_REF] or the equivalent topography, i.e. the geometry plus optionally the roughness coefficient [START_REF] Roux | Use of parameter optimization to estimate a flood wave: Potential applications to remote sensing of rivers[END_REF][START_REF] Castaings | Automatic differentiation: a tool for variational data assimilation and adjoint sensitivity analysis for flood modeling[END_REF][START_REF] Honnorat | Lagrangian data assimilation for river hydraulics simulations Computing and Visualization in[END_REF], mainly from set of sparse trajectories, sparse elevations or a combination of them. The associated dynamic models can consist of simple 1D/2D shallow water optionally coupled together [START_REF] Monnier | On a 2D zoom for 1D shallow-water model: coupling and data assimilation[END_REF] or of more sophisticated models such as the regular ocean modeling system (ROMS) model [START_REF] Powell | 4dvar data assimilation in the intra-americas sea with the regional ocean modeling system (roms)[END_REF] or the Mike11-Mike22 flooding model [START_REF] Madsen | Data assimilation in the MIKE 11 flood forecasting system using kalman filtering in[END_REF][START_REF] Hartnack | Data assimilation in a combined 1D-2D flood model in 'Innovation, advances and implementation of flood forecasting technology[END_REF]. The conceptual bottleneck of DA approaches using noise-free dynamic model is that the behavior of the fluid under study is assumed to be completely explained by tuning the parameters of the model (this is the so-called strong-constraint paradigm of the 4D-Var formulation). However, in practical situations the complete list of parameters involved can be huge: (timedependent) wind stress, (time-dependent) topology, bed roughness, (time-dependent) boundary conditions, ... and for any situation involving parameters not dealt by the data assimilation method, the estimation is likely to fail. This in particular explains that such methods have not been exploited in the purpose of fully characterising the flow under study, with the notable exception of [START_REF] Strub | Inverse estimation of open boundary conditions in tidal channels[END_REF].

Data-assimilation using stochastic dynamic model

An alternative, and to our opinion a more appropriate, track consists in considering that we have only an estimate of the true dynamic model describing the evolution of the flow under study (in the context 4D-Var technique, this is the so-called weakconstraint paradigm). To model this view, one considers that the actual state of the system at each time t results from the propagation of the state at time t -1 through the estimated (known) dynamic plus a stochastic term. By using such a modeling, one is tempted to consider the problem of fully estimating the flow understudy (instead of a set of parameters of the dynamic model). Some authors have followed this track and most of them uses ensemble-based methods [START_REF] Salman | A method for assimilating lagrangian data into a shallow-water-equation ocean model[END_REF][START_REF] Verlaan | Data Assimilation Schemes for Non-Linear Shallow Water Flow Models in[END_REF][START_REF] Madsen | Data assimilation in the MIKE 11 flood forecasting system using kalman filtering in[END_REF][START_REF] Hartnack | Data assimilation in a combined 1D-2D flood model in 'Innovation, advances and implementation of flood forecasting technology[END_REF][START_REF] Vreugdenhil | Application of data assimilation to threedimensional hydrodynamics: The case of the Odra Lagoonn in[END_REF] to solve the problem.

Among them, some authors considered the assimilation of a small set of lagrangian observations, given for example by drifters (Salman et al 2006, Strub, Percelay, Tossavainen andBayen 2009). We think that the application of such observations, limited to particular applications as the acquisition of larges set trajectories, can be in practice a complicated task. The other authors proposed to use water level measurements [START_REF] Verlaan | Data Assimilation Schemes for Non-Linear Shallow Water Flow Models in[END_REF][START_REF] Madsen | Data assimilation in the MIKE 11 flood forecasting system using kalman filtering in[END_REF][START_REF] Hartnack | Data assimilation in a combined 1D-2D flood model in 'Innovation, advances and implementation of flood forecasting technology[END_REF][START_REF] Vreugdenhil | Application of data assimilation to threedimensional hydrodynamics: The case of the Odra Lagoonn in[END_REF]. In these works, the authors consider the assimilation of a very few measurements and are interested on the reliability of the observed components of the flow (velocity or elevation) at some critical points and do not consider the opportunity to characterize the full flow understudy. Two main reasons can explain this situation. In the above mentioned articles and in many hydrological applications, observations consist of sparse local measure of velocity, flow or elevation from which it seems difficult to estimate a reliable complete characterization of the flow all over the space. Furthermore, the use of dense observations necessarily add a significant computational cost and memory usage into the DA procedure. Thanks to the advent of depth-sensors and of Monte Carlo methods for high-dimensional problems, one can reconsider this question.

In this work, we propose a complete system to characterize free-surface shallow flows. This system is composed of a very affordable commercial depth camera and of the particle filter data-assimilation method proposed by [START_REF] Papadakis | Data assimilation with the Weighted Ensemble Kalman Filter[END_REF], dealing with stochastic dynamic model and involved in this study to account for the complexity of real free-surface flows. In section 2, we present the stochastic data assimilation method reconstructing the spatiotemporal geometry and motion fields of free-surface shallow flows from sequences of depth images, only. We propose an enhancement of the particle filter for unknown inflow boundaries by using two observations instead of one. In section 4.1, we discuss the ability of the Kinect depth sensor having the strong advantages to be cheap and of easy use. Then, in section 3 and in section 4, we evaluate the potential of our system by showing its ability to estimate free-surface shallow flows spatiotemporal geometry and velocity fields from synthetic and real data, respectively. Finally, we conclude and give some perspectives in section 5.

The particle filter to assimilate depth images into a shallow flow model

The problem of estimating the consecutive states x 1...t of a physical system from a sequence of partial observations y 1...t is ubiquitous in geosciences. A common approach consists in stating it as a filtering problem i.e. characterizing the probability distributions associated to the state space of the physical system given the past observations and a stochastic dynamic model. More specifically, one states an observation model and a dynamic model respectively defining the densities p Yt (.|x t ) and p Xt (.|x 0 , . . . x t-1 ) that are then combined via Bayes' rule to define the posterior probability p X 0 ,...Xt (.|y 0 , . . . y t ). This last distribution is called the filtering distribution and generally cannot be given in a computable closed-form solution but can be estimated using Monte Carlo methods.

Stating the filtering distribution

The overall model is graphically represented in figure 1 and presented in details in the three next subsections.

The dynamic model

The shallow-water model (often referred to as Saint-Venant model) is practically suited to handle many hydrological situations such as tides, storm surges, river and coastal flows, lake flows, tsunamis and more generally is theoretically supported by the linear wave theory (or airy wave theory) when the ratio water depth under wave length is smaller than 5.10 -2 [START_REF] Vreugdenhil | Application of data assimilation to threedimensional hydrodynamics: The case of the Odra Lagoonn in[END_REF]. Let briefly summarise it. We consider a liquid in a gravitational field in open channel (or free surface) contained in a rectangular cuboid of length L. The free surface level is supposed to be a time and space differentiable function h(x, y, t). The bed topology is flat as sketched in figure 2.1.1. In addition, we assume that the fluid is only subject to gravity, viscous stress and

h 0 u -∆t→0 v -∆t→0 - f x init h ∆t u 0→∆t v 0→∆t - f ... - f h k∆t u (k-1)∆t→k∆t v (k-1)∆t→k∆t - f ? g y 1 ... - f h 2k∆t u (2k-1)∆t→2k∆t v (2k-1)∆t→2k∆t - f ? g y 2 ... ... Figure 1
. Schematic view of the modeling. Each state of the system consists of the elevation and the velocity components for the given time k∆t. The evolution of the system is handled by a function f which is unknown (we simply know a simplified version f ). An observation y t is generated from the states each k time step of the dynamic. The unknown quantities to estimate (i.e. the states x t ) consist in the first row of this scheme while the known variables consist in the second row.

that the liquid is Newtonian, incompressible and of constant viscosity ν. The shallow water approximation consists in considering that = h max /L 1 (calling h max the maximal water level) and that the dimensionless parameters follow the asymptotical dominance as → 0: /Re = O( ), 1/( Re) = O(1) and F r 2 = O(1), with Re = u max h max /ν the Reynolds number and F r = u max / √ gh max the Froude number. Then, performing a dimensional analysis of the Navier-Stokes equations for an incompressible flow according to these dominances, integrating the velocity components u and v along the water depth [0, h], considering that the pressure is hydrostatic and empirically replacing the two bottom stress component ν ∂u/∂z(x, y, 0) and ν ∂v/∂z(x, y, 0) by (gn 2 /h 1/3 )u|(u, v)| and (gn 2 /h 1/3 )v|(u, v)| (from Groves and Groen formula and Manning formula, n being the Manning bed roughness factor, u and v being described in a few lines), the shallow water equations read

∂h ∂t + ∂hu ∂x + ∂hv ∂y = 0 (1a) ∂(hu) ∂t + ∂ ∂x hu 2 + 1 2 gh 2 + ∂(huv) ∂y + gn 2 h 1/3 u|(u, v)| = 0 (1b) ∂(hv) ∂t + ∂(huv) ∂x + ∂ ∂y hv 2 + 1 2 gh 2 + gn 2 h 1/3 v|(u, v)| = 0, ( 1c
)
where u = h(x,y,t) 0 u(x, y, z)dz and v = h(x,y,t) 0 v(x, y, z)dz are the depth integrated horizontal velocities.

On that basis, we consider that the dynamic of the free surface flow under study is completely described at time t by the state variable x t = [u(x, y, t), v(x, y, t), h(x, y, t)], satisfying the discrete time integration of the shallow water system dynamics f :

x t = f (x t-∆t ).
(2)

As previously mentioned, this numerical model is based on physical assumptions and on numerical truncation errors. In other terms, it is attempt to produce only an approximation of the actual dynamic of the flow. We model this aspect by considering the evolution of the state as a stochastic process (in the following we use capital letters to denote random vectors and caligraphic letters to denote matrices):

X t = f (X t-∆t ) + W f t , (3)
with X 0 ∼ N (x init , R 0 ) and W f t ∼ N (0, R t ) (that denotes that the random variable W f t has a normal distribution with mean 0 and covariance matrix R t ). In practice, we model the spatial components of W f t as a stationnary random field whose covariance coefficients evolve as a function of exp(-r 2 /r 2 h ) where r is the distance between the sites of the field and r h is a bandwith parameter [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using monte carlo methods to forecast error statistics[END_REF]. The correlation between the different times and between the three components h, u and v of W f t is chosen null. We call (3) the dynamic model and we note p Xt (.|x t-∆t ) the subsequent conditional probability density function that we call prior. We consider that the discretetime variables X 0 , X ∆t , X 2∆t , . . . X T are connected in a first order Markov chain (i.e. p(x t |x t-∆t , x t-2∆t , . . . x 0 ) = p(x t |x t-∆t )) with an initial distribution p X 0 (.). Each of these variables X t is latent (or unobserved) and fully characterizes the state of the system under study at time t. We call n x the size of each variable X t .

The observation model

We consider that we acquire an observation y t of the geometry of the free surface z = h(x, y, t) using a perspective range camera (often refer to as 2.5D camera). We call g the function such that y t = g(x t ) = g(h(., ., t)). There is no simple formula giving the value of the observation y t for a given elevation function h. A first solution consists in using a ray casting procedure to simulate the operator g. However, such a solution can be very time consuming considering the number of observations that we will have to perform to solve the filtering problem. To simulate the range camera operator, the observation y t can be back-projected into the state grid x t using an appropriate matrix G converting camera frame coordinates to object frame coordinates. This way g(x t ) simply consists in subsampling the elevation component of x t . We tested both methods and we do not noticed any significant difference. For that reason, we used the back-projection approach in order to speed-up the computations.

In the following, the values y t are considered to be rearranged in a vector form (to facilitate the algebraic formulations). As the acquisition process is subject to noise, we model the observation of the state x t as a stochastic process given by

Y t = g(X t ) + W g t , (4)
where W g t ∼ N (0, Q t ). In practice, we consider that there are no spatial correlation, thus Q t is a diagonal matrix. We call (4) the observation model. We consider the set of discrete-time variables Y 1 , . . . , Y t , . . . , Y T conditionally independent provided that X 0 . . . X t , . . . , X T are known. We note p Yt (.|x t ) the subsequent conditional probability density function. We call n y the size of each variable Y t .

The filtering model

The overall model formally reads

     X 0 ∼ N (x init , R 0 ) X t |X t-∆t = x t-∆t ∼ N (f (x t-∆t ), R t ) Y t |X t = x t ∼ N (g(x t ), Q t ).
(5)

Then, given an initial guess of the system at time t = 0 (x init ,R 0 ) and the set of observations y 1,...,T , we want to estimate the filtering laws p X 0...t (.|y 1...t ) (for t ∈ [0, T ]) or one of their estimator such as the root mean square estimator

Ê(X 0...t ) = p(x 0...t |y 1...t )x 0...t dx 0...t . (6)
In the present study this estimator provides the elevation ĥ and the depth integrated horizontal velocities û and v.

Solving the filtering problem : the Weighted Ensemble Kalman Filter

From a theoretical point of view, the evolution of the conditional density p X 0...t (.|y 1...t ) is described by the Fokker-Planck (or forward Kolmogorov) equation [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]). However, its direct numerical evaluation is intractable in real situations. In this context, the representation of the conditional pdf as an ensemble or a set of particles and the integration through time of each of the ensemble elements leads to particularly tractable algorithms. The sequential importance sampling (SIS) filters [START_REF] Doucet | On sequential monte carlo sampling methods for bayesian filtering[END_REF][START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation Radar and Signal Processing[END_REF] and the ensemble Kalman filters (EnKF) [START_REF] Evensen | The Ensemble Kalman Filter: theoretical formulation and practical implementation[END_REF], including its numerous variants, are probably the most used methods. In this work, we followed the recent works of [START_REF] Papadakis | Data assimilation with the Weighted Ensemble Kalman Filter[END_REF] who proposed the Weighted Ensemble Kalman Filter (WEnKF), a sequential importance sampling filter relying on an ensemble Kalman filter. It combines the best of particle filter and EnKF. The algorithm is based on the Ensemble Kalman Filter updates of samples in order to define a proposal distribution for the particle filter that depends on the history of measurement. In the present study we propose to enhanced the WEnKF scheme with the use of two observations y t and y t+1 instead of one. This provides a better coupling of elevation and velocity component estimations when considering unknown inflow boundary conditions.

In the following, we simplify the notations with ∆t = 1 so that e.g. X t+1 = X t+∆t . In section 2.2.1, we recall the foundation of the sequential importance sampling filter. In section 2.2.2, we discuss the use of an ensemble Kalman to build a specific sequential importance sampling filter. In section 2.2.3, we propose to modify the stochastic process underlying the ensemble Kalman filter procedure so that it uses two observations y t and y t+1 instead of one. In section 2.2.4, we discuss a method to build efficiently the Kalman gain, which is the main ingredient of the ensemble Kalman filter.

2.2.1. Bayesian filtering of a non-linear process using sequential importance sampling The sequential importance sampling method stands on the point that we can compute any estimator of p X 0...t (.|y 1...t ) using sequences of samples (x i t ) iteratively drawn from a proposal density function q (optionally parametrized in y 1...T ). We denote this general proposal as q X Q t (.|x 0,...,t-1 ; y 1...T ) (calling X Q t the associated random vector). In particular, we can write the root mean square estimator (RMSE) associated to the filtering distribution as

Ê(X 0...t ) = lim N →∞ 1 N i w t (x i 0...t )x i 0...t 1 N i w t (x i 0...t ) , (7)
where each sample

x i k (k ∈ [1, . . . , t], i ∈ [1, . . . , N ]) is generated from the distribution of density q X Q k (.|x i 0,...,k-1 ; y 1...T
) and where the weights w t are built recursively as

w t (x i 0...t ) = w t-1 (x i 0...t-1 ) p(y t |x i t )p(x i t |x i t-1 ) q(x i t |x i 0..t-1 ; y 1...T ) , (8)
with for t = 0, w 0 (x i 0 ) = 1/N and each sample x i 0 is drawn according to the distribution of density p X 0 (.).

The combination of the equations ( 7) and ( 8) is the keystone of any SIS filters and leads to the following algorithm:

Algo SIS: basic sequential importance sampling algorithm init:

for i = 1 : N w i 1 = 1 N draw x i 0 with density p X 0 for t = 1 : T Sampling: for i = 1 : N draw x i t with density q X Q t (.|x i 0,...,t-1 ; y 1...T ) Weighting: for i = 1 : N w i t = w i t-1 (p(y t |x i t )p(x i t |x i t-1 ))/q(x i t |x i 0:t-1 ; y 1...T ) Normalising: for i = 1 : N ŵi t = w i t / j w j t Ê[X 0...t ] = i ŵi t x i 0...t .
2.2.2. The proposal distribution q In practice, any proposal q X Q t (.|x t-1,...,0 ; y 1...T ) can be used and the different instanciations of the SIS filter mainly differ by this choice. There are two main characteristics to consider when choosing the proposal q X Q t (.|x t-1,...,0 ; y 1...T ). Firstly, it must be possible to sample it a moderate time (sampling q is needed for the SIS algorithm). Secondly, it has to be as close as possible to the target pdf p Xt (.|y 1...t ) to provide a reliable estimation of the target estimator (for a given number of particle N ). A classical choice for q consists of q X Q t (.|x t-1,...,0 ; y 1...T ) = p Xt (.|x t-1 ), resulting in the so-called bootstrap algorithm [START_REF] Gordon | Novel approach to nonlinear/non-Gaussian Bayesian state estimation Radar and Signal Processing[END_REF]. Another classical choice consists of q X Q t (.|x t-1,...,0 ; y 1...T ) = p Xt (.|x t-1 , y t ), resulting in the socalled optimal SIS algorithm [START_REF] Doucet | On sequential monte carlo sampling methods for bayesian filtering[END_REF]. While theoretically supported, these instances of the SIS algorithm are often considered as practically less efficient than the ensemble Kalman filter. However it seems apropriate to mention that unlike the SIS algorithm, the ensemble Kalman filter does not converge, as the number of ensemble members tends to infinity, to the target estimator. In this work, we combine the benefits of the two algorithms by using an ensemble Kalman filter to build the proposal distributions q of a SIS filter.

Interestingly, Le [START_REF] Gland | Large Sample Asymptotics for the Ensemble Kalman Filter[END_REF] showed that the ensemble Kalman filter can be seen as a Monte Carlo simulation of the following stochastic process:

X Q t = X Q f t + K t (y t + W g t -g(X Q f t )), with X Q f t = f (X Q t-1 ) + W f t ( f
standing for forecast) and where K t the Kalman gain associated to samples

X Q f t reads K t = cov(X Q f t , g(X Q f t ))(R t + cov(g(X Q f t ), g(X Q f t ))) -1 .
In practice, the Kalman gain K t is estimated from the set of samples used in the simulation process. As a consequence the simulated random variables X Q t are not independent and identically distributed. However, it has been shown by [START_REF] Gland | Large Sample Asymptotics for the Ensemble Kalman Filter[END_REF] that when N tends to infinity, K t tends to its deterministic limits and the simulated random variables tend to be independent and identically distributed according to the normal distribution:

X Q t |X Q t-1 = x t-1 ; y 1,...,T ∼ (9) N (f (x t-1 ) + K t (y t -g(f (x t-1 ))), (I + K t G)R t (I + K t G) T + K t Q t K T t
), G being the matrix associated to the observation operator g (see Section 2.1.2). This finally permits to define properly a proposal distribution q X Q t (.|x t-1,...,0 ; y 1...T ) and a procedure (the ensemble Kalman filter) to sample it (in the limit of large samples). Then, one can easily restate Algo SIS as:

Algo WEnKF: Sequential importance sampling with Kalman proposal

for i = 1 : N draw x i 0 with density p X 0 w i 0 = 1 N for t = 1 . . . T 1) Draw x i t : 1.1) Forecast step: for i = 1 : N x i f t = f (x i t-1 ) + w i t with w i t generated from N (0, Q t ) 1.2) Analysis step: for i = 1 : N y i t = y t + v i t with v i t generated from N (0, R t ) build K t = cov((x f i t ), (g(x f i t ))) × (cov((g(x f i t )), (g(x f i t ))) + R t ) -1 for i = 1 : N x i t = x i f t + K t (y i t -g(x i f t )) set x i 0...t = [x i 0...t-1 , x i t ] 2) Weight samples: for i = 1 : N w i t = w i t-1 (p(y t |x i t )p(x i t |x i t-1 ))/q(x i t |x i t-1 ; y t ) 3) Compute estimators: Ê[x 1...t ] = i (w i t / j w j t )x i 1...t
A resulting challenge consists in performing the step 2) of this algorithm : evaluating the three densities p Yt (.|x t ), p Xt (.|y t ) and q X Q t (.|x t-1 ; y t ). Whereas p Yt (.|x t ), can be evaluated very efficiently, the two other terms are more difficult to address. Following previous works [START_REF] Hoteit | Particle Kalman filtering for data assimilation in Meteorology and Oceanography[END_REF][START_REF] Papadakis | Data assimilation with the Weighted Ensemble Kalman Filter[END_REF], we chose to drop these two terms. More justifications about this choice can be found in [START_REF] Papadakis | Data assimilation with the Weighted Ensemble Kalman Filter[END_REF]. Further work is needed to evaluate the theoretical impact of this choice.

Using two observations for a better coupling of the state variables

To improve the previously stated algorithm, we slightly modify the proposal distribution q to incorporate two observations y t and y t+1 instead of one. By doing this, we do not simply consider the relationship linking the observation of the current state g(x t ) and the state x t itself to compute the Kalman gain but also introduce the next observation y t+1 . Our motivation is to provide a better coupling between the estimation of the elevation and the velocity components of the state variables by analyzing through the observation y t and y t+1 their combined effects over the time interval [t, t + 1] instead of a simply considering successively each observation independently from each others. In practice, this yields to the following stochastic process:

X Q t = X Q f t + K t,t+1 ([y t , y t+1 ] + [W g t , W g t+1 ] -[g(X Q f t ), g(X Q f f t+1 )]), (10)
where [v 1 , v 2 ] denotes the concatenation of vector v 1 and v 2 and with,

X Q f t = f (X Q t-1 ) + W f t , X Q f f t+1 = f (X Q f t ) + W f f t , with W f f t ∼ N (0, R t )
and where K t,t+1 is the Kalman gain expressing the linear relationship between the observations y t and y t+1 and the state variable X Q f t . It writes:

K t,t+1 = cov(X Q f t , [g(X Q f t ), g(X Q f f t+1 )]) × cov [g(X Q f t ), g(X Q f f t+1 )], [g(X Q f t ), g(X Q f f t+1 )] + R 2ny -1
, where R 2ny is the 2n y × 2n y block diagonal matrix composed of two copies of the matrix R. We expect this process to give more suited samples to approximate the filtering distribution p X 0...T (.|y 1,...T ) and we use it as a generator for the proposal distribution q.

The question of whether this process properly defines a density q is not addressed in this work and is a future research direction. As observed in the validation section, such a choice improved quantitatively the estimations.

2.2.4. Computing the Kalman gain K t To our opinion, the main challenge in the EnKF procedure consists in computing efficiently reliable Kalman gain K t,t+1 (or K t in the original version). Indeed, this matrix is crucial as it determines how to convert the errors from the observation grid into the state grid and ultimately how to compute the correction term associated to x i t regarding the observation y t and the past uncertainty. We assume the covariance matrix R t associated to the observation noise as known and hence focus our interest on how to estimate the covariance matrices cov

X Q f t , [g(X Q f t ), g(X Q f f t+1 )] and cov [g(X Q f t ), g(X Q f f t+1 )], [g(X Q f t ), g(X Q f f t+1 )] (or cov X Q f t , g(X Q f t ) and cov g(X Q f t ), g(X Q f t )
in the basic version). In the original EnKF implementation, these covariance matrices are estimated as the empirical covariance matrices of the samples §. However, the empirical covariance is only a rough estimate of the true covariance when N << n x (this is a case of the so-called large p-small n inference problems). Indeed, the number of degree of freedom (n x (n x + 1)/2) is too high to express correctly the covariance matrices using only N samples. However, X t (and g(X t )) represents a space-structured random vectors and we can infer a structure on their covariance matrix. Adding a structure to or regularizing the covariance matrices will decrease the number of degree of freedom of the statistic. Then, the estimated regularised covariance matrices are likely to better characterize the sample statistic (for a given sample size) and thus to result in a more appropriate Kalman gain and analysis. In this context, a pragmatic solution consists in smoothing the covariance matrix and imposing a complete independence between spatially distant points of the state grid (Houtekamer and Mitchell 2001). Such a localization increases virtually the size of the ensemble. In this work, we choose to use a strategy similar to the one proposed by Houtekamer and Mitchell (2001), that consists in performing the product between the empirical covariance coefficients and a slowly positive definite radial basis function decreasing with point-to-point distance and canceling at a fixed cut-off distance h correl . In practice, we use a [START_REF] Buhmann | Radial Basis Functions[END_REF] function. Then all matrices involved in the analysis step (step 1.2 in Algo SIS-KF) become sparse and the analysis step becomes a sparse linear algebra problem that can be solved efficiently using appropriate optimization algorithms [START_REF] Golub | Matrix Computations Johns Hopkins Studies in the Mathematical Sciences[END_REF].

Numerical validation

In this section, we investigate the performance of the WEnKF method to estimate the dynamic of a free-surface flow from synthetic sequences of depth data.

The robustness of the proposed method to simulated depth data quality is evaluated for two flow configurations: a collapse of a water column as a toy-example and a flow in an suddenly expanding flume as a more realistic flow. Furthermore, we illustrate the interest of using two observations instead of one observation into the correction step (see section 2.2.3). We note oneObs and twoObs the two WEnKF algorithms using in the correction step one observation and two observations, respectively.

The numerical dynamical model used to simulate the two flow configurations involves a finite volume implementation of the shallow-water system of equations (1a,1b,1c) as proposed by [START_REF] Bradford | Finite-volume model for shallow-water flooding of arbitrary topography[END_REF]. Time integration is performed with a second-order Runge-Kutta scheme. A no-slip boundary condition is applied on the walls.

The true states x1...t are calculated with this numerical model running from an initial condition x0 . Based on that true state we build the synthetic depth data sequences taking every fourty time step ∆t of the numerical model integration. Then, to simulate the depth sensor image quality, we add a white noise of standard deviation σ obs and introduce a percentage of outliers p out in the depth images. That provides the observations y 1...t . Note that for the sake of simplicity the observations are created with the same spatial resolution as the numerical simulation grid.

The initial state x init is obtained by deteriorating the initial true state x0 . For that purpose we add a random noise with a covariance having a large range (as for the noise on the dynamical model) such that the initial state perturbation parameter

E init = || x0 -x init || || x0 || , (11)
is equal to a given quantity (typically 0.1), where ||.|| = √ < ., . > is the L 2 -norm acting on the whole flow domain. The following states x 1...t are estimated by assimilating the observations y 1...t given the perturbated initial state x init .

To assess the accuracy of the free-surface flow reconstruction we compare the estimated flows x1...t to the true ones x1...t such as the assimilation errors read

E û, v = || ŵ -w|| u 0 and E ĥ = || ĥ -h|| h 0 , (12)
where x = ( ŵ, ĥ) = ( û, v, ĥ) and x = ( w, h) = ( ũ, ṽ, h) are the estimated (via assimilation) and true state components, respectively, and where u 0 and h 0 are a characteristic velocity and a characteristic length scale, respectively. In addition, we define the free-run errors

E ū,V = ||w -w|| u 0 and E h = ||h -h|| h 0 , (13)
where x = (w, h) = (ū, v, h) are the free-run state components calculated with the shallow-water numerical model from the initial condition x init up to time t. Then, to give an insight of the gain obtained using the observations driving the model instead of running the model alone, we compute the ratios between assimilation and free-run errors for the velocity E û, v/E ū,v and for the elevation E ĥ/E h . The same set of parameters was used to evaluate our method for both synthetic flow configurations. The number of particles N was set to 100. The initial standard deviations of the covariance matrix R 0 were equal to 0.05 h 0 and 0.25 u 0 for the elevation and the velocity components, respectively. Then the standard deviations of R t were set to 0.04 h 0 and 0.06 u 0 for the elevation and the velocity components, respectively. The bandwith parameter r h indicating the spatial decrease of R 0 and R t covariance coefficients was set to 2 h 0 . The covariance matrix Q t was set diagonal (no spatial dependance on noise) with constant standard deviation equal to 0.114 h 0 (this value was estimated for the Kinect depth sensor, see section 4.1). Finally, h correl , the cut-off distance of the radial basis function used for the localization process described in section 2.2.4, was set equal to 0.6 h 0 .

Collapse of a water column

In this simple case, a small circular column of water placed at the center of a square container collapses under gravity and generates a wave. The initial elevation h 0 of the water column above the initial container water level h ∞ was such that h 0 = 0.5 d 0 , where the water column diameter d 0 = 2 cm and with h ∞ = 1.5 d 0 . The width of the container was equal to L = 10 d 0 . The initial velocity in the container was zero. To non-dimensionnalize the results we define the characteristic time scale t 0 = h 0 /g and the characteristic velocity u 0 = √ gh 0 . The collapse of the water column was simulated with the shallow-water numerical model in a computational domain L x × L y = L × L discretized on a square grid of n x × n y = 200 × 200 points and integrated in time with a time step ∆t u 0 /h 0 = 0.006.

Here, we considered the initial state x init as the "perfect" case described above, i.e. with flat and static initial water surfaces (top of the water column and around in the container), whereas the initial true state x0 , considered as the "real" case, was more complex, i.e. with a large scale random extra component providing non smooth initial water surfaces with local non-zero velocities. The observation y 1...t were built from snapshots of the true state elevation component taken every 40∆t u 0 /h 0 , leading to an observations Strouhal number St obs = h 0 /(40∆t u 0 ) = 4.15.

Figure 3 displays the time evolution of the observations, the free-run simulation, the particle filter estimations and the true at three time steps. The results given by the proposed WEnKF method showed a good agreement with the true for both the elevation and the velocities. As expected for a sequential data assimilation scheme, the agreement improved in time. The free-run states, i.e. the simulation starting from x init , rapidly diverged from the truth whereas the data assimilation scheme, starting from an ensemble of simulations around x init rapidly recovered the true states over time. It should be noted here that the proposed data assimilation approach reconstructs the full states of the free-surface flow (i.e. elevations and velocities) based on depth observations alone. This is quite remarkable, given the limited input information. Following those encouraging results, obtained for parameters fixed to median values, it was legitimate to analyse the sensitivity of the method to the data quality and to the initialisation. To this end we investigated varying parameter horizons, ranging from σ obs = 0.015 d 0 to σ obs = 0.12 d 0 , from p out = 0 to p out = 0.35 and from E init = 0.025 to E init = 0.5.

Figure 4 represents data assimilation errors as a function of the non-dimensional level of noise σ obs /h o , the rate of outliers p out and the rate of initial perturbation E init , for twoObs at t/t 0 = 9.51. Whatever the parameter values, the three estimated state components had the same low error level, below 2%. The rate of outliers (given realistic ratio below 35%) did not influence the estimations. On the contrary the errors increased with the level of noise in the observations. For σ obs /h 0 larger than 10% and 5% the errors reached a plateau slightly decreasing for the elevation h and slightly increasing for the two velocity components ū and v, respectively. Note that this changing slope might be correlated with the noise standard deviations of the observation and dynamical models (covariance matrices R t and Q t ) equal to 11% and 6% of h 0 , respectively. This asymptotic behavior, also reflected in figure 5a, indicated the good robustness of the WEnKF scheme. The slope difference seen between the elevation and the velocity component errors could be due to the fact that the elevation was the observed component. Another feature of the observed component was clearly exhibited when plotting the errors as a function of the initial perturbation level E init (see figures 4 and 5b): the elevation component estimations were not sensitive to E init whereas the velocity component estimation errors slightly increased with E init . Note that the noise standard deviation of the observation model As a conclusion of this first simple flow configuration experiment, the proposed particle filter method exhibited very interesting reconstruction and robustness properties. Furthermore, we indicate that we did not observed any significant differences between the methods oneObs and twoObs. This result can be explained by the fact that in the present configuration, the only source of departure from the free-run flow to the true flow came from the unknown initial state. As a result, one did not need a complex coupling between the model and the observation after the first few time integration because by itself the model was able to explain the flow. The two schemes (oneObs and twoObs) will be compared further more extensively for more realistic flow configurations when adding unknown varying inflow conditions and a stochastic forcing term. In the following we fix data quality and initialisation parameters to p out = 0.1, σ obs = 0.03 d 0 and E init = 0.2, respectively.

Flow in an suddenly expanding flume

In this case study, we considered the flow in a suddenly expanding flume as described in figure 6. The flume consisted of a L = 10 cm long, L wide approach flume before a sudden expansion and a L long, 2L wide expanding flume. The inflow water surface elevation and velocities were oscillatory in phase with a frequency of 1 Hz and amplitudes of 1 cm and 0.22 m/s, respectively. The inflow mean elevation and mean velocity were H in = 1 cm and U in = 0.22 m/s, respectively. The corresponding inflow Froude number F r = U in / √ g H in = 0.7 was lower than unity which indicates that the flow was subcritical, i.e. behaved like a fluvial motion. A characteristic of this flow lies in the sudden expansion where the flow separates and generates a two dimensional vortex. In addition we considered two cases for the oscillatory inflow: one with a uniform inlet velocity profile and another with a half bell-shape (squarecosinus) inlet velocity profile leading to a more complex free surface flow dynamic. Both expanding flume flow configurations were simulated with the shallow-water numerical model in a computational domain L x × L y = 2L × 2L discretized on a square grid of n x × n y = 200 × 200 points and with a time step ∆t u 0 /L = 0.006, where u 0 = √ g H in . The true inflow states xin,0...t were varying periodically in elevation and in velocity. For the assimilation methods, the variations in time of the inflow conditions were considered as unknown and in practice were fixed to their mean values H in and U in . We stress the fact that to show the ability of the stochastic model (3) to handle unexpected situations, the inflow components were not considered as an unknown vector parameter of an augmented model system like in [START_REF] Gronskis | Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation[END_REF], but just as state variable inlet values. Moreover we considered the case, named stochastic forcing, where we add a random large scale extra component to the state vectors when generating true states x0...t before each time integration. Hence, each true states included the stochastic forcing transported by the dynamical model during one time step ∆t and more. This was done to mimic the complexity of real free surface flows which cannot be simulated with the deterministic shallow water model (1a,1b,1c). In contrast with data assimilation methods considering "perfect" models, the stochastic data-assimilation scheme (5) carried out in this study better modeled such complex situations. The noisy observations y 1...t were built from the true state elevation component taken every 40∆t u 0 /L, thus leading to an observations Strouhal number St obs = L/(40∆t u 0 ) = 4.17.

Results showing the particle filter estimations for an homogeneous inlet velocity profile are given in figure 7 together with the true states. The free surface elevation and motion were reconstructed well for all the time steps considered, but were in less agreement near the inlet area especially when the surface dynamic was more complex due to the stochastic forcing. As reflected in figure 8 this behaviour was emphasized with more complex half-bell shape inlet velocity profile. The errors of the estimations performed in the configuration with stochastic forcing and uniform inlet velocity profile are displayed in figure 9. In the inlet region the highest errors reached 20% and 43% for the elevation and for the motion, respectively. However, the use of two observations (twoObs) instead of only one (oneObs) improved significantly the agreement for both the elevation and the motion. This demonstrated the capability of the proposed particle filter scheme to deal with complex free surface flow configurations, and as it was expected the use of more observations was necessary to reconstruct better the flow. Table 1 summarizes the influence of the inlet velocity profile shape, of the stochastic forcing and of the number of observations considered, on the estimations accuracy. These results also indicate the gain obtained using the observations driving the model instead of running the model alone. The ratios between assimilation and free-run errors for the velocity E Û , V /E U,V and for the elevation E ĥ/E h increased with the complexity of the flow and decreased when using two observations instead of one. This illustrates the importance of coupling the data and the dynamical model to reconstruct the overall flow state (elevation and motion) with a single partial observation (elevation at the inflow location at a given time).

Experimental demonstration

We now characterize the ability of the proposed WEnKF scheme, to reconstruct free surface real flows geometry and motion, given only Kinect based measurements of the elevation.

Kinect depth sensor

Capturing the geometry of a 3D object is of great interest for many applications. Different strategies (e.g. stereovision, shape from deformations) handling different technologies (e.g. binocular system, structured light) have lead to a large variety of sensors having different characteristics. We report the reader to the work of [START_REF] Herbort | An introduction to image-based 3D surface reconstruction and a survey of photometric stereo methods 3[END_REF] for an introduction to depth sensors. Note that recently stereoscopic particle images velocimetry approaches have been proposed to estimate simultaneously morphology and velocity of moving free surfaces [START_REF] Turney | A method for three-dimensional interfacial particle image velocimetry (3D-IPIV) of an airwater interface[END_REF][START_REF] Gomit | Free surface measurement by stereo-refraction[END_REF], Chatellier et al 2013). In this work, we choose to use the Kinect sensor to estimate the geometry of a free surface flow for the first time [START_REF] Combès | Free-surface flows from Kinect: feasability and limits in[END_REF]. More recently [START_REF] Mankoff | The Kinect: a low-cost, high-resolution, short-range 3D camera Earth Surf[END_REF] introduced the Kinect sensor to the earth science community. This depth sensor has the advantage to be cheap (about 150 e), to provide high-level programming interfaces and not to need any specific calibration. These three characteristics make it a good candidate for practical use. It is partly composed of a 320 × 240 RGB sensor at 30 Hz, a 320 × 240 infrared sensor at 30 Hz and an infrared pattern projector. Range images are obtained from the so-called light-coding technique: the projection of the infrared pattern on the object under study is captured by the infrared sensor and the analysis of this projection is used to recover the geometry of the object. This approach is close to the optical profilometry technique proposed by [START_REF] Cobelli | Global measurement of water waves by Fourier transform profilometry[END_REF] for the measurement of water waves and in practice provides range images of 640 × 480 pixels at 30 Hz. To the best of our knowledge, the Kinect sensor is not well-documented and the main available information are shared by users from the web. Note that in the paper of [START_REF] Mankoff | The Kinect: a low-cost, high-resolution, short-range 3D camera Earth Surf[END_REF] hardware and software are described and a code for data processing is provided.

In this work, we do not focus on the technology of the Kinect and of its subsequent limitations and facilities but use it, as best as we can, in a "black-box" way. For that purpose, we investigated its accuracy to estimate range images of solid and liquid smooth surfaces. Liquid's light diffusivity was enhanced by the addition of white dye for the use of the Kinect sensor. For a sake of simplicity, our study was restricted to surfaces located between 680 and 780 mm from the device. However, note that the devise is able to estimate depth data up to 13 meters leading to an observable surface of about 10 m × 14 m with a magnification of about 2 cm near the optical axis and 20 cm at the periphery of the observable window. A more complete description of our experiments can be found in [START_REF] Combès | Free-surface flows from Kinect: feasability and limits in[END_REF]. The main conclusion of our studies were the following: On solid surface, the Kinect sensor displayed a measurement uncertainty of 0.9 mm for both flat and sinus-like surfaces; The sensor captured successfully sinus-like varying elevations with spatial periods smaller than 20 mm and amplitudes smaller than 2 mm; Measurement errors coming from observations of solid and liquid surfaces were comparable when the attenuation coefficient of the liquid was larger than 113 m -1 . As an illustration, figure 10 shows a 3D temporal reconstruction of a water wave moving in a rectangular flat bottom tank, captured by the Kinect.

Wave in a rectangular flat bottom tank

The real experiments carried out in this study consisted in observing the free surface of a fluid contained in a rectangular flat bottom tank of size Lx × Ly = 250 mm × 100 mm as illustrated in figure 10. More specifically we observed the evolution of a unidirectional wave generated by an initial free surface height difference h 0 = 1 cm. In the following the characteristic velocity u 0 is considered as an approximation of wave phase velocity √ g h 0 . The wave propagation was simulated with the shallow-water numerical model in a computational domain L x × L y discretized on a square grid of n x × n y = 222 × 88 points and with a time step ∆t u 0 /L x = 0.0042.

The Kinect depth observations were characterized by a high level of noise and exhibited large regions of missing data on the boundaries due to light reflections on the tanks wall. The observed free surface behaved roughly as an unidirectional wave along the x axis. The initial state x init was considered with a flat surface and a null velocity although the real state was a moving wave. The observations y 1...t were assimilated every 10∆t u 0 /L, leading to an observations Strouhal number St obs = L/(10∆t u 0 ) 24, that was rather high. We tested ensemble size from 20 to 200 and in practice when N > 100 we do no more observe any significant reduction the estimation error. As a result, we use N = 100. Assimilating each observation took about 10 minutes on a standard personal computer.

Figure 11 shows estimated free-surface elevation colormap and motion vector fields for a water wave propagating in a rectangular flat bottom tank. This time sequence in a plan view exhibits one back and forth of the wave. The reconstructed height sequence is a compromise between noisy Kinect observations and an imperfect dynamical model, exhibiting the spatiotemporal complexity of the free surface motion. The predicted free surface motion fields indicated the highest velocities on the wave crest and recirculation regions on either side upstream and downstream of the wave. Providing a physically consistent initial state for such a real flow was not straightforward: we could use the measured elevation however the associated motion fields was difficult to model. Interestingly, although we initiated the assimilation process with a flat surface and a null velocity, as time proceeds the particle filter rapidly provided accurate estimations. The whole free surface geometry dynamics was remarkably reconstructed from noisy measures of the elevation only.

Conclusion & perspectives

Most of the time there are only temporal sequences of sparse data coming from lagrangian or eulerian measurements of the flow elevation to characterize a free surface flow dynamics. Although sparse and noisy, these observations when properly coupled with a shallow water model contain enough information to reconstruct the whole free surface flow dynamics.

This paper has presented the first free surface flow reconstruction using Kinect based depth measurements only. The Kinect sensor was able to capture observations of wave-like surfaces with wavelengths and amplitudes sufficiently small to support applications such as flow monitoring or medium to large scale flows characterization. The sequential data assimilation algorithm, based on a particle filter stochastic approach, has been validated on two numerical cases and a real experiment in laboratory. Results have shown remarkable reconstruction of both elevation and velocities of the free surface flows. In addition, the influence of noise and outliers in the input depth data has been assessed, indicating that the method has exhibited good robustness properties up to large deteriorations of the data quality. Another feature of the proposed WEnKF scheme has been the quasi non-sensitivity to the initial state perturbation levels. When starting from flat surface and null velocity the technique has rapidly and accurately recovered the dynamics of complex free surface flows. The enhancement of the particle filter with two observations instead of classically using only one has improved significantly the agreement for both the elevation and the motion. Finally, in contrast with data assimilation methods considering "perfect" models, the stochastic data-assimilation scheme carried out in this study has yielded remarkable estimations, given the complexity of the flow configurations and the limitations of the dynamical model.

Note that all the presented results characterized acquisitions within an area near the image center and from a short range. Future works will consist in investigating the data assimilation technique over larger ranges (we recall that the Kinect is able to estimate depth data up to 13 meters leading to an observable surface of about 10 m × 14 m with a magnification of about 2 cm near the optical axis and 20 cm at the periphery of the observable window) and under more practical conditions such as discharge estimation. Future works will also consist in estimating parameters of the model such as the bed roughness in addition to the model noise term in order to reduce the departure of our model to the actual dynamic of the fluid and ultimately to estimate more accurately the states of the system.

Figure 2 .

 2 Figure 2. Free surface flow configuration.

Figure 3 .

 3 Figure3. Comparisons of elevation colormaps in millimeters and of free-surface motion vector fields, for the water column collapse. From top to bottom: observations y, free-run states x, assimilated states x, true states x. From left to right: t/t 0 = 3.17, t/t 0 = 6.34 and t/t 0 = 9.51. Perturbation parameters on the observations and on the initial state were set to σ obs /h 0 = 0.16, p out = 0.1 and E init = 0.1.

Figure 4 .

 4 Figure4. Results for the collapse of a water column, for twoObs at t/t 0 = 9.51: Left, elevation component errors E ĥ; Right, velocity components errors E û, v as defined in(12). From top to bottom, errors as a function of: the non-dimensional level of noise σ obs /h o with p out = 0.1 and E init = 0.1; the rate of outliers p out with σ obs /h 0 = 0.1 and E init = 0.1; the rate of initial perturbation E init with p out = 0.1 and σ obs /h 0 = 0.1.

Figure 10 .

 10 Figure 10. 3D water surface geometry estimated with a Kinect sensor: left, picture of the experimental setting when observing a moving wave in a water tank; right, 3D reconstruction of the wave moving in the tank from Kinect acquisitions at t = 0.15, 0.3, 0.45 and 0.6 seconds. All spatial scales are in millimeters. From Combès et al (2011).

Figure 11 .

 11 Figure 11. Elevation colormaps in millimeters and free-surface motion vector fields for a wave in a rectangular flat bottom tank. From left to right and top to bottom, snapshots of the estimation displayed every 20 time step ∆tu 0 /L x .

Table 1 .

 1 Ratios between assimilation and free-run errors for the suddenly expanding flume flow at time t u 0 /L = 1.98. E û, v /E ū,v E ĥ/E h E û, v /E ū,v included in the true states x0...t and therefore in the observations y 1...t .

	oneObs	twoObs

a

§ This choice leads to the algebrical reformulation exploited by the reduced rank, and the EnTKF formulation.