Assimilating SAR-derived water level data into a hydraulic model: a case study

Satellite-based active microwave sensors not only provide synoptic overviews of flooded areas, but also offer an effective way to estimate spatially distributed river water levels. If rapidly produced and processed, these data can be used for updating hydraulic models in near real-time. The usefulness of such approaches with real event data sets provided by currently existing sensors has yet to be demonstrated. In this case study, a Particle Filter-based assimilation scheme is used to integrate ERS-2 SAR and ENVISAT ASAR-derived water level data into a one-dimensional (1-D) hydraulic model of the Alzette River. Two variants of the Particle Filter assimilation scheme are proposed with a global and local particle weighting procedure. The first option finds the best water stage line across all cross sections, while the second option finds the best solution at individual cross sections. The variant that is to be preferred depends on the level of confidence that is attributed to the observations or to the model. The results show that the Particle Filter-based assimilation of remote sensing-derived water elevation data provides a significant reduction in the uncertainty at the analysis step. Moreover, it is shown that the periodical updating of hydraulic models through the proposed assimilation scheme leads to an improvement of model predictions over several time steps. However, the performance of the assimilation depends on the skill of the hydraulic model and the quality of the observation data. Correspondence to: L. Giustarini (giustari@lippmann.lu)


Introduction
Due to its all weather and day and night capability, Synthetic Aperture Radar (SAR) is regarded as the most promising technology to monitor floods from space. Since the launch of the ENVISAT mission in 2002 and more recently the successful launches of the high-resolution COSMO Skymed, TerraSAR-X and Radarsat-2 missions in 2007, considerable progress has been made with respect to SAR-based flood delineation algorithms (e.g. Zwenzner and Voigt, 2009;Martinis et al., 2009;Mason et al., 2010;Matgen et al., 2010Matgen et al., , 2011. These methods were specifically developed for rapid, repeatable and reliable flood mapping. Remote sensing data have become more frequent and rapidly available and accuracies of SAR-derived flood detection have improved due to higher spatial resolutions and enhanced image processing algorithms. There is a growing pressure on the scientific community to find new ways to use the increased volume and accuracy of remote sensing data in order to improve near realtime flood monitoring and prediction applications (Di Baldassarre et al., 2009).
The retrieval of water level data by merging remote sensing-derived shorelines with a digital elevation model ("indirect measuring technique") can be viewed as a way to add value to remote sensing data for hydrological applications (e.g. Hostache et al., 2009;Raclot, 2006;. Direct measuring techniques such as those from the proposed swath altimetry "Surface Water and Ocean Topography" (SWOT) mission  represent a potential enhancement of the indirect measuring techniques as they enable the systematic acquisition of elevation data Published by Copernicus Publications on behalf of the European Geosciences Union.
L. Giustarini et al.: Assimilating SAR-derived water level data of inland water surfaces with an observation uncertainty of 50 cm for a 50 m spatial resolution (Lee et al., 2010). Both techniques enable the monitoring of changes in water volume in ways that are not possible using hydrometric station data.
However, as space-borne sensors provide instantaneous snap-shots of an area of the Earth's surface, there is a need to combine remote sensing data sets with hydrologic-hydraulic prediction models to generate time-lapses of flooded surfaces. Sequential data assimilation methods can be used to integrate time-continuous model state forecasts (e.g. soil moisture, surface water storage) with remote sensing observations as they become available. The quantification of uncertainty for all observations is a pre-requisite to any meaningful data assimilation study. To date, only a few studies have investigated the uncertainty description of remote sensing derived water stages. According to , geo-location accuracy of the flood extent and parameter uncertainty in flood delineation algorithms are the most significant sources of uncertainty in a high resolution ENVISAT ASAR-based flood mapping application. In a water level retrieval process, which consists in merging the remote sensing-derived shorelines with a digital elevation model, these errors add up to the errors that are inherent in the topography data. The uncertainty assessment approach of  results in cross-section specific cumulative distribution functions (cdfs) of water elevation estimates. In the case study of the Alzette 2003 flood, the indirect water stage measuring technique yields cdfs that indicate non-normal distributions and skewness of the SARbased water level estimates for many cross-sections. Uncertainty of stage over the entire river reach was on the order of 2 m.
In situ measurements are routinely assimilated for hydrologic-hydraulic modelling applications (e.g. Neal et al., 2007;Madsen and Skotner, 2005;Pauwels and De Lannoy, 2009). However, only a few studies have attempted to assimilate remote sensing-derived water stage data into hydraulic models. Matgen et al. (2007) used a direct insertion method that forced water stage data simulated by a hydraulic model to fall within the confidence interval of ENVISAT and European Remote Sensing satellite (ERS-2) SAR-derived water stages. They showed that the insertion of remote sensing derived water levels increased the accuracy of modelled water levels. However, this version of a direct insertion method is not an optimal sequential assimilation method and appears as a useful approach only if uncertainties associated with observation data are much smaller than simulation uncertainties and distribution functions of observations are unknown.
Different variants of the Kalman filter present dynamic methods for merging uncertain simulation data with uncertain observations. Andreadis et al. (2007) and Biancamaria et al. (2010) successfully applied an ensemble Kalman filter (EnKF) to assimilate synthetic water level measurements from the proposed SWOT mission with simulations from the LISFLOOD-FP 2-D hydraulic model (Bates and De Roo, 2000). Durand et al. (2008) assimilated virtual SWOTderived water level observations into a hydraulic model of the Amazon River to improve the estimates of bathymetric depths by 84 % compared to the model runs without assimilation. Neal et al. (2009) used the EnKF with a real event SAR image of the flood extent. The ensemble uncertainty was estimated by image histogram thresholding with different backscattering values and repeatedly shifting the resulting flood boundaries in space in order to approximate geolocation errors. The measurement error covariance was defined from the perturbations of this ensemble of water level estimates around the mean. Neal et al. (2009) only considered the measurement members with the smallest interquantile range over the ensemble at any river section. This was suggested as a quality control step prior to assimilation that was needed because some locations produced biased data (e.g. shorelines next to steep slopes and tall vegetation). They also argue that, because of the spatial coverage offered by remote sensing, it is not necessary to use all measurements. Following this approach, they showed that it is possible to significantly reduce discharge and water level uncertainty of a hydraulic model by using ENVISAT Advanced SAR-derived water stage estimates.
Since the Gaussian error assumption may not be satisfied for most remote sensing observations of water stage, Matgen et al. (2010) proposed an assimilation scheme based on the Particle Filter (PF) as a possibility to relax the Gaussian assumption in the EnKF while preserving its advantages. Their experiments showed that the PF is able to correct water depths from a corrupted 1-D hydrodynamic model by assimilating synthetic observations that are realistic in terms of accuracy for remote sensing-derived water levels. In this case, the PF leads to a significant increase of the accuracy and a reduction of the model forecast uncertainty. Matgen et al. (2010) further state that problems related to a spatially and temporally variable non-Gaussian distribution of water level observations still need to be solved.
The objective of this paper is thus to examine the usefulness of currently available satellite data to update a hydraulic model in near real time, through a PF-based assimilation scheme. The specific objectives are: (1) to adapt the PF assimilation scheme in order to deal with non-Gaussian distributions of remote sensing derived (RSD) water levels; (2) to deal with model structural errors and parameter uncertainties, proposing two variants of the PF; (3) to assess the usefulness of SAR data with respect to in situ hydrometric station data.

Study area and available data
The area of interest is situated in the Grand Duchy of Luxembourg (Fig. 1 The hydrologic model is applied to a basin area of 356 km 2 draining to a stream gauge located in Pfaffenthal. This produces the upstream boundary condition for the hydrodynamic model, which simulates the 19 km reach of the Alzette River between hydrometric stations at Pfaffenthal and Mersch (Montanari et al., 2009). The river reach is described by 144 ground-surveyed and evenly spaced (∼130 m) channel cross sections.
The investigated event was a flood recorded in January 2003. Hydrometric data were recorded every 15 min at six stream gauges interspersed throughout the 19 km reach (Pfaffenthal, Walferdange, Steinsel, Hunsdorf, Lintgen and Mersch). Moreover, information about the maximum water level reached along the river during the flood was available, as measured by means of a theodolite (altimetric accuracy around ±2 cm) at distributed points across the floodplain. The availability of in situ data not only allows evaluating the results of the assimilation of remote sensing data, but also helps to contrast the use of space-based and in situ based water level monitoring in a data assimilation exercise. The comparison of results provides insights on the advantages and limitations of each data set.
This paper makes use of the measured precipitation rate to drive the hydrologic model: hourly rainfall data observed in Livange between 1 and 7 January 2003 are available. Contrary to Neal et al. (2009), who used predictions of convective and stratiform precipitation and evaporation, in this case the forcings of the hydrologic model can be considered as the best available representation of the rain field, potentially leading to a more accurate basin response. The set-up of this case study can be viewed as a realistic representation of an operational application of the proposed methodology. Two subsequent SAR images, acquired at two distinct stages of the 2003 flood event have been used in this study: one was acquired by the ERS-2 satellite during the rising limb of the flood wave; the second image by the ENVISAT satellite just after the flood peak. The two images are shown in Fig. 2 as well as the six stream gauges distributed along the river. A LiDAR DEM of the floodplain at a spatial resolution of 2 m and a vertical accuracy of ±15 cm was fused with remote sensing derived flood boundaries to retrieve the water stages . RSD observations of the water stages in the river were retrieved from the two available images, which have a spatial resolution of 12.5 m: in other words, at each cross section more or less independent water stages were observed. L. Giustarini et al.: Assimilating SAR-derived water level data Uncertainty assessment for RSD water stages Different approaches can be used to estimate spatially distributed water levels and their associated uncertainties from a sequence of wet/dry flood edges extracted from radar imagery and fused with a DEM.  proposed a procedure to estimate uncertainty associated with RSD water stage data using a Monte Carlo-based statistical analysis.  introduced a slightly different approach based on a more integrated uncertainty assessment and based on the analysis of the confidence that can be given to the SAR derived shorelines. Both approaches take into account different sources of uncertainty (i.e. parameters of image segmentation algorithm, co-registration of geo-information layers, accuracy of digital elevation model) that affect the retrieval of water elevation data from remote sensing imagery. First, flood extension limits with their respective geo-location uncertainty are derived from a SAR image using a radiometric thresholding-based procedure. Next, the part of the SAR derived shoreline having a high probability of being erroneous (close to building and trees) due to the incapability of the SAR sensor to detect water in the corresponding areas are removed from the analysis. This provide pertinent shorelines that will be used for the water level estimation. Finally, the ensemble of relevant flood boundaries is superimposed on a DEM in order to estimate water levels. The method takes into account the uncertainty stemming from the underlying DEM and ultimately provides empirical distribution functions of water level data from space at every river cross section (Fig. 3). In assimilation studies, this approach thus potentially allows exploiting the full empirical distribution of observed water levels. However, the resulting uncertainty can be very high (Fig. 3). Moreover the distribution functions often exhibit bias and skewness, especially in the vicinity of steep embankments. All of these factors render the use of the empirical distribution functions in data assimilation studies problematic.
In order to reduce the estimation uncertainty, all water level estimates were hydraulically constrained. The procedure was first introduced by Raclot (2006) and consisted in applying hydraulic rules, governing overland flow in a floodplain, on water level intervals derived from aerial photographs. The hypothesis for applying this procedure is that the uncertainty on the estimated water level is known accurately enough to assume that the real water level is included inside the water level estimate intervals. Considering the effort that was made to take account of all sources of uncertainty and remove errors impacting the extraction of water levels from remote sensing imagery, it is reasonable to assume that the "true" water level is included inside each interval. This hypothesis is supported by the fact that all ground-surveyed measurements of water elevation are included in the above-mentioned intervals . Up-/downstream relationships between water level estimates were first defined depending on the location of the water stage estimates within the floodplain. Knowing that the water level decreases from upstream to downstream, an algorithm imposed the following two constraints on the water stage estimate intervals: (1) the upper bounds of the water stage intervals have to decrease from upstream to downstream and (2) the lower bounds of the water stage intervals have to increase from downstream to upstream. This algorithm allowed a significant reduction of the mean water level estimation intervals, as shown in the panels on the bottom right in Fig. 3. As a result, water level information was available as cross-section specific values of the possible local water levels.
3 Simulation design Figure 4 shows the setup of the data assimilation experiments using event data. The methodology consists of assimilating remote sensing-derived water stage observations into an ensemble of 1-D hydraulic model integrations for a number of cross sections. The upstream boundary conditions (flow hydrographs) are produced using an ensemble of semi-distributed hydrologic model forecasts with perturbed parameter sets, initial conditions and precipitation data.

Coupled hydrologic-hydraulic model
The semi-distributed hydrologic model is loosely coupled to a 1-D hydraulic model: the discharge hydrograph computed by the hydrologic model is used as upstream boundary condition to drive the hydraulic simulation, but the hydraulic model does not feed back into the hydrologic model.
The rainfall-runoff transformation is carried out using the Community Land Model version 2.0 (CLM 2.0) (Dai et al., 2003), a global land surface model built over the 356 km 2 drainage area of the Alzette River extending upstream from the gauging station at Pfaffenthal. To generate meaningful ensembles of model predictions, we followed the procedure of De Lannoy et al. (2007). Model parameters, forcings and initial conditions of the hydrologic model were perturbed in such a way that the ensemble mean differs from the observation by a value that is equal to the time average of the L. Giustarini et al.: Assimilating SAR-derived water level data ensemble spread (De Lannoy et al., 2006). More details on the ensemble generation and the verification measures that were used to monitor the ensemble spread can be found in Matgen et al. (2010). The land surface model was initialized a month before the analyzed flood event to allow spin up and balancing of the state variables in each ensemble member.
The assumption is made that model uncertainties derive only from the upstream boundary condition, which means that uncertainties in hydraulic model structure, parameters, geometry and lateral inflow are not accounted for. Moreover, it is worth mentioning that errors in the inflow peak timing are not taken into account. Therefore, in order to represent the hydrodynamic model uncertainty, an ensemble of 64 upstream boundary conditions (i.e. discharge) was generated with CLM 2.0, adopting the same methodology as discussed in Matgen et al. (2010). The hydrographs are shown in Fig. 5 together with the 2 time steps of satellite overpasses for the 2003 flood event. As opposed to the synthetic experiment by Matgen et al. (2010), no artificial bias was added here to the output of the hydrologic model. However, we noticed that during the receding limb of the hydrograph, the ensemble did not bracket the observed discharge, indicating a tendency of the hydrologic model to underestimate the inflow during that period. The ensemble of the hydraulic model realizations or particles (see box on top-left of Fig. 4) has been produced by integrating the hydrodynamic model with all the members of the discharge ensemble generated by the hydrologic model for the analysis period 1 January 15:00-7 January 23:00, 2003 (GMT+1).
The hydraulic model is implemented over a 19 km reach of the Alzette River between the gauging stations of Pfaffenthal and Mersch. Since flow direction in this area is mainly parallel to the channel, the 2-D flow field that is typically related to riverbank overtopping can be accurately approximated by a 1-D representation and thus, the Hydrologic Engineering Center River Analysis System -HEC-RAS (HEC-RAS 4.0, 2008) -was set-up for 1-D river flow computation.

Data assimilation algorithm
The data assimilation technique applied in this study is a sequential Particle Filter (PF), an ensemble filtering method that has its origin in Bayesian estimation. Unlike the widely used EnKF (e.g. Burgers et al., 1998;Evensen, 1994), which simplifies the recursive estimation by assuming a Gaussian distribution for both the model and the observation error structure, the PF relaxes the need for restrictive assumptions regarding the shape of the probability density functions and can easily manage the propagation of a non-Gaussian distribution through nonlinear hydrologic and hydrodynamic models (Moradkhani, 2008). In the PF the assumption of Gaussianity is relaxed and the fundamental idea is to represent the required posterior density by a set of properly weighted samples (Smith et al., 2008), named particles, and to compute the estimate based on these samples and weights.

Sequential Importance Sampling (SIS)
The simplest implementation of the PF is based on the Sequential Importance Sampling (SIS) method (see box in the middle of Fig. 4). Each particle consists of a possible value of the state. The filtering density is approximated by a discrete distribution, whose support is the set of particles. The probability mass assigned to each particle is proportional to that particle's weight, which, in turn, is proportional to the likelihood of the observation at the assimilation time step (Fearnhead, 2002). In this case study, a particle represents the water surface line resulting from one hydrodynamic model run at the assimilation time step t = k, and the number of state variables corresponds to the number of cross sections. The particles are sampled directly from the state-space to represent the posterior probability, and a weight is computed for each particle according to the information contained in the RSD water level observations. In this case study, the SIS algorithm was implemented using two different distribution functions that are characteristic for the data sets at hand. A local weight, w i,j , is assigned to any state variable j for any particle i (the index k is left out to avoid confusion but the weights are recomputed at any time-step of observation acquisition). Note that the weighting procedure can be adapted to any kind of distribution function.
Section 2.1 describes an approach for retrieving for each cross section an interval of possible water stage values. In this case, as only the maximum and minimum water level values are available for each cross section, we assimilate the observation data assuming a uniform distribution, whose boundaries are the derived maximum and the minimum water stage estimates. Therefore, the likelihood or weight, w i,j , of the water level, x i,j , simulated by particle i at cross section j for the time-step of the observation acquisition is simply computed as: where x is the state vector of the state variables (simulated water stages at any cross section j for any particle i),a and b are the lower and upper observed endpoints or minimum and maximum, respectively, for any cross section j . All assimilation experiments were also carried out using hydrometric station data as it allowed contrasting results obtained with relatively uncertain but densely distributed satellite-derived data with those obtained with accurate but poorly distributed in situ measurements. When assimilating data from these stations, a Gaussian distribution was used, assuming the recorded water level to be the mean of a normal distribution whose shape is also defined by a pre-defined value of standard deviation (Matgen et al., 2010). Note that water stage estimates obtained from space-borne swath altimeters were assumed to be normally distributed as well (Andreadis et al., 2007). One weight, w i,j , for any water level, x i,j , simulated by particle i at cross section j for the time-step of the observation acquisition is therefore computed as: where x is the state vector of the state variables (simulated water stages at any cross secrtion j for any particle i), µ is the observation vector and σ is the standard deviation associated to the observations. Regardless of the way the weights of the individual particles are computed (e.g. assuming a normal or uniform distribution of residuals), the matrix of weights contains all local weights, w i,j , that are obtained for any model run or particle i at all the N o observed cross sections j . Subsequently, an overall likelihood of the water level globally simulated along the river reach for any particle or model run is computed by applying the joint probability theory for independent variables: where N o is the number of observations. The resulting global weight is then normalized.
N p is the number of particles or water surface lines.
The probability obtained with the global weights at the previous steps allows computing an expectation of the updated water stage as follows for the assimilation time-step k:

Sequential Importance Resampling (SIR)
Because of the stochastic behaviour of the system, the particles tend toward dispersion and there is a risk that many of them obtain negligible weight during the analysis. This behaviour, called degeneracy, is defined as the tendency to converge to a single point estimate (Moradkhani et al., 2005).
To avoid unnecessarily high computational resources, a resampling step is carried out to eliminate samples with low weight and to replicate samples with high weight. In other words, the HEC-RAS model is re-initialized at the time-step of the observation acquisition, k, replicating the water lines with higher weight (see box on the bottom of Fig. 4). The resampled particles have the same weight, until the next assimilation step. The most common resampling scheme is the Sequential Importance Resampling (SIR) developed by Gordon et al. (1993). The authors refer to Moradkhani et al. (2005) and Weerts et al. (2006) for more detailed explanations of the SIR and its use in hydrologic applications. It should be mentioned that resampling is independent of the proposal distribution (Fearnhead, 2002) but it is also possible to implement a PF without resampling. In any case, the SIR algorithm also suffers from particle degeneracy. Smith et al. (2008) showed that the resampling step only reduces the degeneracy of the particles. Moreover, a different problem may arise, known as sample impoverishment, causing particles with high weight to be selected many times, which leads to a loss of diversity in the sample. In fact, due to the discrete approximation of the filtering density, inaccuracies accumulate over many time steps and the result is often a clustering of particles in small areas of the state-space (Fearnhead, 2002).
The experimental set-up of this case study avoided the problem of sample impoverishment through a loose coupling of the hydrologic and hydraulic model components. Only the water levels were resampled, while the spread in the input data was maintained. Therefore, even in the extreme case where a single particle was retained and replicated, the spread in the discharge hydrographs ensured that shortly after the assimilation a sufficient spread of the state variables was obtained. As in Matgen et al. (2010), at the assimilation time-step k, the estimate of the upstream water level x 1 exp (Eq. 5) was used to compute the corresponding estimate of the discharge E Q 1 k , using the HEC-RAS internal rating curve. A simple algorithm was then applied for the updating of the forcings (i.e. discharge Q) until the next assimilation time-step (see box on bottom right of Fig. 4): where Q 1 i is the upstream discharge of particle i and Q 1 k is the average, computed during the analysis step k, of the hydrograph ensemble. The forcing update was applied until the next assimilation time-step, based on the assumption that relative model errors remain constant and that correcting the inflows by the same relative error term at subsequent time steps will improve the accuracy of the model prediction.

Synthetic experiment vs. real-event case study
The set-up of the case study was very similar to the one presented for the synthetic experiment in Matgen et al. (2010), but there were some significant differences that need to be highlighted.

First assumption: observations
In the synthetic experiment, one of the basic assumptions was that the observed and assimilated (synthetic) satellite water level data were unbiased and normally distributed. In a real case study, the Gaussian assumption may not be satisfied for at least some of the remote sensing-derived water stage observations. Examples of cross-section specific pdfs of RSD data retrieved with the procedure proposed by  are shown in Fig. 6. The data obtained from the first satellite overpass (ERS-2 on 2 January 2003 at 11.00 GMT+1) are given for four representative cross sections. As it can be seen, the data at each section exhibit a different pdf shape. Obviously, the normal distribution is not a suitable candidate distribution for representing the pdf. In the same panels, the maximum and minimum water levels deriving from the hydraulic coherence concept ) are also shown. It can be seen from these results that a significant reduction of the water level estimation intervals was obtained (see also Fig. 3). Here we assume that within each interval the RSD water levels are uniformly distributed.
Furthermore, Fig. 6 highlights bias and skewness in the RSD data retrieved using the first image. The RSD water levels at the first satellite overpass (ERS-2 on 2 January 2003 at 11.00 GMT+1) for the gauged cross section of Lintgen (named 115 in the hydraulic model) were not centered on the in situ water level measurements. For this cross section and the considered time step, the RSD water levels showed a tendency to overestimate the actual water level. After applying the hydraulic coherence concept, the plausible interval was significantly narrowed and, most importantly, included the ground "truth". Therefore, the data assimilation algorithm was run assuming the uniform pdf to represent the statistical distribution of RSD water levels.

Second assumption: model
A further characteristic of the synthetic experiment by Matgen et al. (2010) was that the hydraulic model was correct in it structure, parameter set and initial or analysis conditions. Therefore, the differences between observations and models only derived from inaccuracies in the input data (i.e. hydrographs at the upstream boundary). This means that, for a given forcing (Q), the model generally performs equally well (or poorly) at all cross sections along the river. In other words, due to the fact that the same model is used both to generate the artificial satellite observations and to assimilate them, a model run that is good at a given cross section performs well at all the other cross sections.
In a real case study, model structure errors (e.g. 1-D flow approximation, errors in geometry) and parameter uncertainties (e.g. Manning's roughness values), cause local bias that need to be taken into consideration. We expect such models to have a less uniform behavior along the river reach. This raises difficulties in the selection of a good model run, as it might happen that one model performs globally well over the whole river reach but at the same time has a poor performance at a local level (i.e. at some cross sections).
In the present case study, a first assimilation test was carried out with the same Manning's values as in Montanari et al. (2009): one value for the channel and one for the floodplain. However, by using spatially-distributed stagedischarge measurements it is possible to further reduce potential errors that may originate from a too simplistic parameterization of the model. By introducing additional model parameters, we expect the resulting model to have potentially better behaviour at a local scale.
For the Alzette River and for the considered flood event, previous studies have shown that the floodplain does not play a significant role in the flood hydraulics Montanari et al., 2009). Therefore, for the floodplain a Manning's coefficient equal to 0.184 s m −1/3 was assumed (Montanari et al., 2009). The Manning friction coefficients for the channel was calibrated by means of discharge measurements carried out in the period 2001-2009 along the river reach on the hydrometric stations of Pfaffenthal, Steinsel, Hunsdorf and Lintgen. The geometry of the river is described by the 144 channel cross sections that were surveyed in 2001. It is important to note that the SAR-observed flooding event occurred in January 2003. All cross sections with available simultaneous measurements of water level and discharge were analyzed by comparing the cross section of the hydraulic model with those observed during each discharge measurement campaign. As there is no evidence that would indicate any significant changes in riverbed geometry, we assume the river geometry to be temporally stable.
Measurements of both water levels and discharge taken at the same time allowed calibrating Manning's roughness values for the four cited cross sections in the main channel. Following the random sampling of 4 values for the channel roughness from a range of plausible values (i.e. 0.030-0.060 s m −1/3 ), parameters of the cross sections that are located between the 4 gauging stations were estimated through linear interpolation. As upstream boundary discharge, we used the January 2003 discharge data that was estimated from recorded water levels and a calibrated rating curve. The model was evaluated by comparing the observed rating curves at the 4 cross sections with the internal rating curves of HEC-RAS. The selected model set is the one that minimizes the difference between simulated and observed water levels for a range of observed and simulated discharge values.
The calibration reproduced the high discharge values reasonably well, with Manning's roughness coefficients set equal to 0.042, 0.044, 0.053 and 0.039 s m −1/3 , for Pfaffenthal, Steinsel, Hunsdorf and Lintgen, respectively. The plots in Fig. 7 display the calibrated rating curves for the cross sections with available measurements. For the cross sections in between, the friction coefficient was linearly interpolated. It has to be noted that the distance between Pfaffenthal and Steinsel is nearly 8 km, which is rather significant with respect to the total river reach length (19 km). Therefore we expect model errors to be higher in the upper part of the river reach. In particular, the cross section of Walferdange did not have measurement data and its friction value was deduced from the one in Steinsel, the nearest calibrated cross section.
The observed January 2003 discharge hydrograph at the upstream boundary condition was used as input for the calibrated model to assess its capability in reproducing observed spatio-temporal fluctuations of water surface elevation. This assessment was made by comparing the observed and the simulated hydrographs at all the gauged cross sections. The Nash-Sutcliffe efficiency was computed obtaining an average value on the gauged cross sections of 0.84, indicating a satisfactory result in terms of model performance.

Results and discussion
The hydrodynamic model, built from the 144 surveyed cross sections, was used to simulate water levels along the river reach, considering as inputs the ensemble of 64 hydrographs generated by the CLM 2.0. At each satellite overpass, RSD water level estimates were assimilated into the coupled hydrologic-hydraulic model following the procedure outlined in Sect. 3.2. The results are compared against in situ observed station data. It is important to note that in situ data were assimilated only at the time steps of the satellite overpasses. This analysis can be considered as a benchmark test that enables contrasting the performances obtained when assimilating, respectively, very precise but poorly distributed ground-surveyed information and spatially distributed but highly uncertain satellite data.

Analysis step
When a satellite observation becomes available, weights are computed for all the simulations at any cross section. A global weight is computed for every particle according to the procedure outlined in Sect. 3.2.1. Figure 8 shows the histograms of computed water stages before and after resampling the particles during the assimilation procedure. The histograms shown here correspond to the four intermediate gauging stations in Walferdange, Steinsel, Hunsdorf and Lintgen. The two groups of four panels on top refer to the first satellite overpass (i.e. ERS-2 image), while those on the bottom refer to the second assimilation time step (i.e. ENVISAT image). On the left, the results of the assimilation corresponding to the application of the uniform pdf with the RSD observed data are displayed, while on the right the reported outcome corresponds to the application of the normal pdf with the in situ data. The in situ measurements at six hydrometric stations located on the river reach also serve as validation datasets for assessing the performance of the analysis. The performance of the assimilation is evaluated through the mean error value, the change in distance between the mean of the a priori histogram and the truth compared to the distance between the mean of the a posteriori histogram and the truth. The standard deviation of the histogram both a priori and a posteriori is computed as a measure of the reduction of uncertainty in the water level at the assimilation time step.
Considering the use of the uniform pdf (panels on the left side in Fig. 8), the results obtained via the assimilation of the intervals, defined by the maxima and minima of the retrieved water stages, show a significant reduction of the spread in the a posteriori distribution of the simulated water stages. The reduction in uncertainty is evident for the first time step and becomes even more significant for the second. Moreover, at the time step of the ERS-2 image acquisition, at all the investigated cross sections the a posteriori distribution of water level estimates encompasses the truth. However, the spread reduction is most significant for cross sections in Hunsdorf and Lintgen, located in the downstream part of the river, where most of the flooding occurred (i.e. where most observations of water stage are available). As a result of the analysis there is a decrease of the mean error value, changing from −0.07 to 0.02 m for Hunsdorf cross section and from −0.14 to −0.06 m for Lintgen. The decrease in terms of standard deviation (changing from 0.29 to 0.13 m and from 0.33 to 0.14 m for the two cross sections, respectively) further outlines the positive effect that the assimilation procedure has at a local level. However, for the sections in Walferdange and Steinsel, both located in the upstream part of the river reach, we observe a tendency to overestimate the recorded stage data post assimilation. In both sections the mean error value increases, from 0.22 to 0.37 m in Walferdange and from 0.17 to 0.34 m in Steinsel. The fact that the standard Table 1. Mean error (mev) and standand deviation (std) values computed for the two satellite acquisitions before and after the assimilation analysis step, for the four cross sections with available ground observation measurements, considering the assimilation of the RSD data with the uniform pdf and the use of the normal pdf with the hydrometric station data. deviation is reduced merely illustrates that uncertainty, as it is defined here, only expresses model uncertainty and not truly probabilistic prediction limits of the assimilation procedure. This limitation of the "global weighting procedure" is confirmed by the results obtained at the second assimilation time step. Overall, the assimilation of data retrieved from the ENVISAT ASAR image leads to the selection of only two particles; one is kept for the simulation and the other one is replicated 62 times. This means that only two simulations provide water surface lines that are included in the RSD intervals at all cross sections. As a matter of fact, the standard deviation values are negligible after this time step, while the mean error values are reduced to less than 10 cm for all cross sections, except Walferdange, where an overestimation remains apparent. On all other sections the a posteriori distribution includes the truth. It has to be observed that for Walferdange it was not possible to calibrate the roughness value of the hydraulic model, due to the unavailability of discharge measurements. Its Manning's coefficient was interpolated considering the values of the upstream and downstream cross sections. The poor quality of the model results at this cross section could thus be explained with a badly calibrated model. Moreover, as flooding only occurred on the downstream part of the river reach, there are no RSD observations available in the upper part of the river (i.e. upstream of the Steinsel cross section). These limitations partly explain the difficulty the global weighting procedure has to select models that perform well along the entire river reach. This result shows that one of the main assumption of the synthetic experiment, namely that input data is the only source of error in hydraulic modelling, can no longer be maintained. Table 1 summarizes the results in terms of mean error and standard deviation values for the two assimilation time steps and the two distribution functions.
Before applying and testing a procedure based on local weighting, an experiment was carried out with precise in situ measurements of water level, recorded at six cross sections along the river at the time of the two satellite overpasses. This test represents a circular way to operate, due to the fact that the same data set is used for assimilation and validation. However, the rationale behind this test is to distinguish model errors from observation errors. A good model is expected to provide results that are centred on the truth at every single cross section. In this experiment the normal pdf was used, assuming a standard deviation equal to 0.1 m to represent the uncertainty of the measurements. As can be seen from the panels on the right in Fig. 8, for the ERS-2 satellite overpass, the resampled particles show a good reduction of the spread (standard deviation ranging from 0.04 m to 0.10 m) and always encompass the truth. However, the assimilation of in situ measurements at the time step of the ENVISAT image acquisition highlights a contradictory behaviour. The reduction of the spread is very significant and only a limited number of particles are selected at any cross section. The results are surprisingly similar to those results obtained with the less accurate RSD observations. However, only in Steinsel did the assimilation methodology predict the truth, while for all other cross-sections the resampling of the particles either leads to a slight over-or underestimation. This indicates that no model run performs equally well along the entire river reach. We assume that this is due to the fact that important sources of errors are not sufficiently well represented by the ensemble of model runs. Table 1 summarises the information in terms of performance and standard deviation values for the two satellite acquisitions and the two distribution functions. From these results we conclude that the proposed PF-based filtering approach is an efficient tool to assimilate observations described by characteristic distribution functions. However, the observed over-and underestimations are an indication that local inconsistencies persist in the calibrated model. This can lead to a sub-optimal functioning of the PF or indeed the rejection of all models: when one particle represents the water level over the whole reach as state vector and the weights are computed according to Eq. (3), systematic model errors at a local scale heavily impact the weight that is attributed to individual particles. Rather than selecting and replicating the best particles overall, the application of the joint probability theory for independent particles penalizes particles that have low weight at some locations. As a result, the global weighting procedure favours compromise solutions that provide acceptable results at all model cross-sections. Therefore, it is important to bear in mind that the PF has been initially designed for removing noise and not systematic errors. A pre-requisite for the application of the PF is thus to reduce systematic errors in the model prior to any data assimilation experiment.

Forecasting step
Following the analysis step, the model is propagated in time: to do so, the hydrodynamic model is first re-initialized with updated water stages and then run with the updated upstream inflow data until a new observation becomes available. Figure 9 shows the stage hydrographs corresponding to the cross-section at Lintgen, when using the uniform pdf to assimilate the RSD observations. The performance of the forecast is illustrated considering the RMSE of the ensemble mean water stage, h ass , with respect to the recorded water stage h truth : (7) which is computed over different time windows, starting from the first time step after the satellite overpass, t ass +dt, and stopping at t − (t ass +dt)+1, where dt is the time step of the simulation, in order to evaluate the usefulness of the assimilation as a function of the number of time steps following the analysis step.
As it can be seen from Fig. 9, after the ERS-2 satellite overpass, when the analysis step efficiently drags the simulated water level towards the observed water level, the RMSE is first close to zero before gradually increasing at subsequent time steps due to the predominant effect of the inflow condition. Although the relative error term was correctly inferred from satellite observations, it becomes obvious from Fig. 9 that the proposed inflow correction model (Eq. 6) underpredicts simulation errors in the time window between the two satellite acquisitions. However, the forecasts with filter are consistently better than the open loop predictions for the first time steps. After the second acquisition the RMSE approaches zero and the error terms in Eq. (6) leads to correct Fig. 9. Stage hydrographs at the cross section Lintgen with the 2 assimilation time steps (bottom panel) considering the RSD water level intervals assimilated through a uniform distribution: the forecasting performance is illustrated with the RMSE evolution in time (top panel). The cyan line represents the RMSE before assimilation and the black line displays the RMSE after assimilation. predictions for more than 5 h after the assimilation, as the error remains constant for some time steps. Nevertheless, later on the application of a constant error prediction term for the inflow leads to wrong predictions, with an underestimation of the receding limb of the hydrograph. Hence, the analysis step is of fundamental importance in order to carry out an efficient inflow correction. Errors in the analysis propagate through the inflow correction model and this can significantly decrease model performance at later time steps.

Analysis step
As an alternative to the global weighting procedure, an approach based on local weighting has been developed and tested. In this procedure, each cross section has its own particle set, with as state vector the water levels at the cross section itself (i.e. the state vector is a scalar). Each particle has its own weight, as required for the PF. This procedure is different from global weighting, where one particle has as state vector the water level over the whole reach. Once again, the uniform distribution was used in the assimilation algorithm. We expect a local weighting based method to yield better results when model structural and parameterization errors are present. However, this method will be rather sensitive to local observational errors. The results are given in Fig. 10. The histograms were obtained at four representative cross sections.
Using the uniform distribution together with the local weighting procedure, the reduction of uncertainty is less evident than with the global weighting method. At each cross section all particles that provide water stage estimates included in the RSD intervals are retained and equally weighted. As a result, the standard deviation values are generally higher than before. While the global weighting procedure applied to the ENVISAT-derived data led to the selection of only one particle, the application of a local weighting procedure causes many simulations to be retained at each cross section. For instance, in Lintgen after the second assimilation step the standard deviation is reduced to 0.00 using a global weight but only to 0.24 using local weights. All the a posteriori distributions include the truth, even if in Walferdange a tendency to overestimate the truth persists. With respect to the mean error values, the two weighting approaches give comparable results at the first assimilation step. For the second satellite overpass an improvement can be observed for the upper part of the river, whereas for Hunsdorf and Lintgen, both located in the lower part of the study area, the tendency to slightly underestimate the truth is further enhanced.
Finally, the assimilation was carried out using in situ water level measurements at six cross sections at the time of the two satellite overpasses (Fig. 10). Here we apply the local weighting procedure to compute the a posteriori histograms. As expected, this is the easiest setup for the assimilation algorithm to recognise the truth. For both time steps and in all the cross sections, the resampled particles encompass the truth. This experiment shows that if observations are of high quality, the local weighting procedure yields very satisfactory results. Adopting a local weighting procedure significantly reduced model uncertainty, as demonstrated by the means of the global weighting experiment. However, as we suspect other sources of error than inflows (e.g. spatially varying friction parameters, intermediate inflows and/or errors in geometry) to be responsible for contradicting results obtained in sub-reaches of the model domain, it has to be expected that these improvements cannot be maintained over many time steps. It is therefore recommended to use the results of the analysis to find the reasons for regionally conflicting results and to make use of such a diagnosis for improving the model in a more persistent way than through a mere reinitialization and inflow correction.
Considering the local weighting variant of the PF, Table 2 summarises the results in terms of mean error and standard deviation for the two time satellite acquisitions and for the two distribution functions. Figure 11 illustrates the performance of the forecasts considering the local weighting procedure to assimilate the ground measurements recorded at six hydrometric stations. It shows that when the analysis step gives good results (i.e. when the error term at the time of the observation is correctly estimated), the short-term forecasts with assimilation are improved. However, the limitation of a possible overcorrection of inflow persists as the inflow-corrected mid-term to longterm forecasts seem to have less skill than the open loop predictions. A possible explanation is that the two satellite observations were acquired during the hydrograph's rising limb when model errors are known to be only weakly correlated in time. This is due to the fact that during the rising limb, errors are difficult to predict as precipitation errors continuously add to model parameter and model structural errors. This conclusion is in line with the findings of Matgen et al. (2010) who stated that because of the underlying Fig. 10. Idem Fig. 8, but applying the local weighting approach for every location along the river (lw: local weight). assumption of constant relative errors, the inflow correction model is reliable only during flow recession periods.

Conclusions
Our case study illustrates advantages and drawbacks related to the application, in a quasi-operational context, of a PFbased assimilation of RSD water levels into a hydraulic model. Two variants of the PF, based respectively on local weighting and global weighting procedures, are proposed. In the global weighting procedure, a single particle contains water levels at all cross sections as state vector. Hence, the likelihood for each particle is derived from its ability to correctly predict water levels along the entire river reach. The local weighting procedure attributes a separate particle set to each cross section (i.e. a single particle has the water level from one cross section as state vector) and thus associates likelihoods to each particle according to its ability to correctly predict water stage at a given cross section. The experiment concludes with the following findings.
1. Matgen et al. (2010) demonstrated through a series of synthetic experiments with unbiased model forecasts and observations that a PF-based assimilation scheme enables the sequentially updating of flood forecasting models. The filter helps to correct for errors in the forcings and guides the recovery of the correct water depth over a modelled river reach. In our real-event case study, even according to the best-case scenario when precise in situ measurements are assimilated into a hydraulic model, difficulties arise from the fact that model accuracy varies in space. This makes it difficult for a global weighting procedure to identify a model run that performs equally well at all cross sections. In fact, input errors are not the only source of model uncertainty. Parameter uncertainty and geometry errors, as well as intermediate inflow errors, lead to locally biased model results. Therefore, an assimilation scheme that is based on a local weighting procedure seems to be the preferred solution when dealing with a model that cannot be well calibrated and when observations with a very low uncertainty are to be assimilated.
2. Before any assimilation of data, the set-up and calibration of the hydraulic model are of paramount importance. It is important to bear in mind that the Particle Fig. 11. Idem Fig. 9, but for the local weighting procedure and the assimilation of the ground measurements.
Filter (and also the Kalman Filter) is a method designed to filter noise, not systematic errors. Our results show that this is particularly important when in situ measurements are assimilated, because there is a significant risk that the performance of a model at a local level is not truly representative for its behaviour at a regional level. In this case, the assimilation can potentially lead to a deterioration of model performance.
3. The quality of the observation data is the second factor that largely determines the effectiveness of the filter. In instrumented basins with well-calibrated models, the hydrodynamic model uncertainty appears low compared to currently available remote sensing observation uncertainty. Nevertheless, our experiment shows that forecast improvements are achievable with currently existing SAR data. However, there is a need to take into account the possibility of bias in the observation data. A PF that is based on a local weighting procedure is the preferred solution when assimilating unbiased and/or very precise observations as it helps to identify the true water surface line at the time of data acquisition. In ungauged basins where RSD water levels are the only available data source, the PF with a global weighting procedure is to be recommended. Certainly, methodologies for retrieving water levels from remote sensing observations need to be improved. The availability of VHR SAR satellites and the global DEMs with increased accuracy can be used to further reduce uncertainty and bias of such data sets. The hydraulic coherence concept that was applied in this study is another step forward. We show that observational uncertainty can be significantly reduced by using hydraulic rules governing overland flow in a floodplain to correct unrealistic water levels.
through data assimilation procedures. Moreover, the application of the proposed assimilation scheme to other case studies will lead to a better understanding of the scaling issue linking the length of the river reach to the model forecast performance.