Benefits and limitations of data assimilation for discharge forecasting using an event-based rainfall-runoff model

Mediterranean catchments in southern France are threatened by potentially devastating fast floods which are difficult to anticipate. In order to improve the skill of rainfallrunoff models in predicting such flash floods, hydrologists use data assimilation techniques to provide real-time updates of the model using observational data. This approach seeks to reduce the uncertainties present in different components of the hydrological model (forcing, parameters or state variables) in order to minimize the error in simulated discharges. This article presents a data assimilation procedure, the best linear unbiased estimator (BLUE), used with the goal of improving the peak discharge predictions generated by an event-based hydrological model Soil Conservation Service lag and route (SCS-LR). For a given prediction date, selected model inputs are corrected by assimilating discharge data observed at the basin outlet. This study is conducted on the Lez Mediterranean basin in southern France. The key objectives of this article are (i) to select the parameter(s) which allow for the most efficient and reliable correction of the simulated discharges, (ii) to demonstrate the impact of the correction of the initial condition upon simulated discharges, and (iii) to identify and understand conditions in which this technique fails to improve the forecast skill. The correction of the initial moisture deficit of the soil reservoir proves to be the most efficient control parameter for adjusting the peak discharge. Using data assimilation, this correction leads to an average of 12 % improvement in the flood peak magnitude forecast in 75 % of cases. The investigation of the other 25 % of cases points out a number of precautions for the appropriate use of this data assimilation procedure.


Introduction
Mediterranean catchments in southern France are exposed to intense rain events that may result in devastating flash floods such as the 2002 event in the Gard Department or that of the Aude Department in 1999 (Gaume et al., 2009). In order to better forecast these events, hydrologists may look to rainfall-runoff models. These may be either physically based, like MARINE  or CASC2D (O'Donnel, 2002), conceptual, reservoir-based (Kobold and Kay, 2004;Bailly-Comte et al., 2012;Tramblay et al., 2010;Fleury et al., 2007) or "black box" type (Toukourou et al., 2009;Piotrowski et al., 2006). All these models are subject to a number of uncertainties tied either to their structure, which implies a simplification of actual physical processes, or to the dataset used for external forcing or the calibration of parameters.
A data assimilation procedure can be applied in order take into account such uncertainties. Within the framework of inverse problem theory (Tarantola, 2004), this approach is intended to combine two types of information -one coming from observations and the other from a numerical model -in order to propose an improved estimation of the system state (Bouttier and Courtier, 1999). First introduced in meteorology (Daley, 1991) and oceanography (Ghil and Malanotte-Rizzoli, 1991), data assimilation techniques have become more widespread in the geosciences, following the increased availability of remote sensing data (McLaughlin, 2002;Reichle, 2008), as well as in hydrogeology (De Marsily et al., 1999) or hydraulics Ricci et al., 2011). In hydrology, assimilation techniques are typically employed in flood forecasting and are referred to as updating techniques. These techniques seek to correct the components of a hydrological model (forcing, parameters, state variables, or discharges) in order to improve the quality of discharge predictions (Refsgaard, 1997) as observational data become available. The most widely used techniques remain those relying on autoregressive (AR) models in order to correct the simulated discharge output by the model (Yang and Michel, 2000;Xiong and O'Connor, 2002). This type of correction is based on the structure of the error between observed and simulated discharges but does not take into account the source of uncertainty. Certain techniques focus on correcting input variables (such as precipitations) (Khal and Nachtnebel, 2008), which often constitute the primary source of uncertainty in operational forecasting. However, this approach has only been infrequently used (Moore et al., 2005). Other techniques emphasize correcting the hydrological model's state variables either "empirically" (Weisse et al., 2003) or more formally via a best linear unbiased estimator (BLUE) (Thirel et al., 2010) or a Kalman filter algorithm (Da Ros and Borga, 1997;Aubert et al., 2003). In  and Maradhkani and Hsu (2005) it was shown that a dual state-parameter estimation using either an ensemble Kalman filter or a particle filter properly accounted for uncertainties in model inputs, outputs and parameters. Variational techniques have also been used to correct hydrological model parameters (Yang and Michel, 2000;Bessière et al., 2007).
The study presented here is aimed at improving, through use of data assimilation, the discharge forecast produced by an event-based hydrological model. Starting with the model state for a given date, the procedure will seek to correct model parameters using existing observations, in this case discharges recorded at the basin outlet. This correction will be calculated using a simplified Kalman filter algorithm, known as best linear unbiased estimator (BLUE). The study is conducted on a Mediterranean basin located in the South of France. The objectives of this article are threefold: (i) to select the parameter(s) that offer the most efficient correction of the peak discharge to be forecasted; (ii) to demonstrate the impact of the correction of the initial moisture deficit of the soil reservoir on simulated discharges at the catchment outlet using data assimilation; and (iii) to identify and understand situations when the assimilation procedure leads to a degradation of the simulated discharge. Section 2 presents the hydrological model as well as the assimilation method and hydrological dataset used to simulate the flood events on which the procedure is tested. Section 3 offers a sensitivity study of model parameters, along with the results of the data assimilation procedure implemented to correct the parameters, which have the greatest impact on model results. Section 4 provides an analysis of cases where the assimilation method was not effective and proposes a number of precautions to be applied when using this technique.

The hydrological model
The simulation of discharges for this study is performed with an event-based, distributed, parsimonious rainfall-runoff model. It has been implemented on a set of regular and independent grid cells (water does not flow between grid cells, but rather is transported directly to the outlet). This model combines a Soil Conservation Service derived (SCSderived) runoff function with a "lag-and-route" routing function. SCS-LR model was presented in Coustau et al. (2012), and will be briefly summarized as follows: i. The catchment is divided into a mesh of regular and independent grid cells.
ii. The rainfall is interpolated in each cell at each time step.
iii. The total runoff is determined by using two reservoirs (cf. Fig. 1 in Coustau et al., 2012) -both empty at the beginning of the event. The cumulated rainfall reservoir has an infinite capacity and determines the direct runoff using a runoff coefficient derived from the SCS runoff model, as used by Gaume et al. (2004): Table 1. Selected flood events -QH p : discharge peak with an hourly time step (m 3 s −1 ); rt: response time (the delay between rainfall peak and discharge peak) (h); P : average cumulative rainfall (mm) calculated according to the Thiessen method; S cal : the value (mm) of the initial condition after calibration; Nash cal : the Nash value after calibration; Hu2 ini the value of the Hu2 index used to initialize the initial condition of the model. The values shown in bold correspond to simulations run with rainfall radar data. where l m is the distance between the cell m and the outlet, V and K 0 are 2 routing parameters respectively acting as the routing speed and the diffusion coefficient.
v. The sum of the elementary hydrographs from all the cells of the catchment, at all the time steps of the event provides the complete hydrograph of the flood.
The complete hydrological model features a total of 5 parameters: S, w, ds, V and K 0 , which remain uniform over all catchment grid cells. This model was implemented using the ATHYS modelling platform (www.athys-soft.org).

Available data
This study was conducted on the Lez Mediterranean catchment in southern France; this catchment covers an area of 114 km 2 and is located upstream of the city of Montpellier. Floods primarily occur during autumn and winter and are caused by intense rainfall events (sometimes reaching several hundreds of mm within a 24-h period), characterized by short response times (on the order of 2 to 5 h) with high peak discharge (up to 480 m 3 s −1 in September 2005). For this study 12 flood events are selected between 1994 and 2008. All of these events display a peak discharge above 90 m 3 s −1 , which corresponds to a return period of 2 yr or more as displayed in Table 1.
It has been demonstrated that, in the Lez Basin, radar rainfall at the beginning of autumn is of higher quality than at the end of autumn or during winter. Since the basin is located at a distance of approximately 60 km from the radar station, the limited vertical extension of clouds coupled with the low elevation of the 0 • C isotherm during winter results in a poor quality of the measurements .
Radar rainfall data were used to simulate the events taking place on October 2001, September 2003, September 2005and October 2008. Ground rainfall data were used for all other events. Ground precipitation data are provided with an hourly time step and recorded using 4 rain gauges -one located on the topographic basin and the other three within a radius of 5 to 10 km outside the basin (cf. Fig. 3 in Coustau et al., 2012). Ground precipitation data were interpolated according to the Thiessen polygon method. Between 1994 and 2000, only the Prades rain gauge on the topographic basin determined the rainfall. Since 2000 the Saint-Martin and Montpellier rain gauges have allowed for a better representation of the rainfall, notably in the west and the south of the catchment. The Mauguio rain gauge is used when a data gap occurs for one of the other three rain gauges.

Model calibration and initialization
The five parameters (S, w, ds, V and K 0 ) of the model were calibrated using 21 events (including the 12 events of Table 1) observed between 1994 and 2008 . These parameters are first calibrated separately for each event. Next, a global calibration procedure, which consists on averaging the value of each parameter over the entire set of events, was applied for all parameters, except for S, which represents the initial water deficit at the beginning of the event. Thus, S varies from one event to another. The optimal S value (denoted S cal in Table 1) is determined for each event by maximizing the Nash criterion where Q sim (t) are the simulated discharges, Q obs (t) the observed discharges and Q obs the average of observed  discharges. For the 12 events selected for data assimilation experiments, the Nash values (denoted Nash cal in Table 1), obtained with S cal and the batch-calibrated parameters (ds = 0.28 day −1 ; w = 101 mm; V = 1.3 m s −1 ; and K 0 = 0.3), range between 0.60 and 0.94 with a median value of 0.86. Thus, using an optimal initialization, the model is shown to satisfactorily reproduce the 12 selected flood events. However, the S cal value cannot be determined while the event is going on and, in an operational situation, the S value has to be set with external predictors of the initial wetness state of the basin. The S cal values are correlated with various predictors including the root-layer humidity output by the Safran-Isba-Modcou (SIM) distributed hydrometeorological model (Habets et al., 2008;Quintana Seguí et al., 2009) (denoted Hu2 ini in Table 1) and the piezometric level, both taken at the beginning of the event. Linear regressions established between these predictors and S cal are used to estimate a value of S, denoted S reg , and allow for the initialization of the event-based model in an operational framework . For example, S reg is calculated from Hu2 ini according to the following equation: S reg = −8.84 Hu2 ini + 732 (R 2 = 0.69). The values of S reg used for the model initialization are presented in Table 2. They correspond to the background value S b used by the data assimilation method.

The data assimilation method
The data assimilation method used in this study is the best linear unbiased estimator algorithm. In this case, the most sensitive parameters of the model are corrected using the discharge observations at the catchment outlet. The correction is calculated over a single time window and not cycled. Because of the non-linear relationship between the simulated discharges and the most sensitive parameters of the hydrological model, an outer loop was necessary to overcome this difficulty. The variational assimilation method (Bouttier and Courtier, 1999) consists of minimizing the cost function, J (x) in order to obtain the optimal values of the variables stored in a control vector, x: The control vector, length n, contains the set of n variables to be optimized. These variables may correspond to the forcing, state variables, or model parameters. x b is the background vector of length n, which contains the a priori values of the control vector variables. y o is a vector of length p, which contains the p observations to be assimilated. The cost function J (x) quantifies both the distance between the control vector x and the background x b , and the distance between the control vector x projected in the observation space and the observations y o , weighted respectively by the error covariance matrices B and R. B, with a size of n × n, and R, of size p × p, are both symmetric positive definite matrices containing the covariances of the error in the background x b and in the observations y o , respectively. The errors in x b and y o are assumed to be independent, Gaussian and unbiased.
H is the observation operator: it maps the control space onto the observation space and is usually non-linear. As a consequence, the cost function J (x) is not quadratic and hence difficult to minimize. To overcome this difficulty, J (x) is locally approximated using an incremental approach (Courtier et al., 1994) which leads to a quadratic cost function, J inc l (δx) that is easier to minimize: where δx = x − x b and H l is derived from the linearization of the observation operator about a reference vector, x l . H l is approximated here by the forward finite difference for a small perturbation δx = (..., δx i , ...) such that its i-th column reads This function J inc l (δx) is minimum when its gradient is zero: In order to solve ∇J inc l (δx a l ) = 0 for δx a l , an iterative process that can be applied using a minimizer or the solution can be derived analytically as in the best linear unbiased estimator (BLUE) algorithm. Since the size of the error covariance matrices B and R is small in the present case, the direct inversion of the BLUE algorithm (Gelb, 1974;Talagrand, 1997;Bouttier and Courtier, 1999, Sect. 4) was chosen to find the minimum of the cost function, J inc l (δx). This minimum can be calculated as where K l is the gain matrix (Eq. 10) (Bouttier and Courtier, 1999), d l is the innovation vector (Eq. 14) representing the difference between the observation vector y o and the linear approximation of the control vector in the observation space The error covariance matrix of x a l is expressed as When calculating the BLUE analysis, x a l , the incremental approximation J inc l (δx) is minimized rather than the true cost function, J (x). In order to approximate the minimum of the true cost function, as done so with the 3D-VAR approach, an outer loop is applied to the BLUE algorithm. This loop iteratively updates the linearization point value by setting the background equal to the analysis of the previous iteration (x l+1 = x a l ). During the first iteration of this outer loop, the linearization step is performed around the background x 1 = x b , whereas during subsequent iterations the linearization occurs around the analysis from the previous iteration. This approach accounts for some of the non-linearities in the observation operator, H .

Sensitivity study to model parameters
The assimilation procedure uses discharge observations at the catchment outlet in order to correct hydrological model parameters after the calibration and initialization steps. A sensitivity study of simulated discharges to various model parameters was carried out; this study serves to identify those parameters which most heavily influence the flood discharge arrival time and magnitude as well as those most suitable for correction with the data assimilation algorithm.
The sensitivity study was conducted on the October 2001 event, which is representative of major flood events occurring in the study catchment (Fig. 1). The observed hydrograph (blue curve) indicates a single flood peak of approximately 240 m 3 s −1 . The simplicity of the hydrograph shape makes it possible to accurately identify the influence of each model parameter on discharges at the outlet. A "free run" (black curve) was integrated with the calibrated parameter values w, ds, V and K, in addition to S reg = 160 mm, derived from the linear regression established with the Hu2 indicator. This simulation shows an underestimation of the peak flow, which can be explained by various errors: uncertainty in the rainfall estimation, error in the initialization or error in the model structure. Despite an underestimation of the flood peak, this simulation provides a satisfactory depiction of discharges observed at the outlet with a Nash criterion value of 0.88.
Starting from the free run, five "perturbed" simulations were conducted by perturbing each of the five parameters, one at a time, by +10 % with respect to its calibrated value. In order to quantify the influence of these parameters on discharges at the catchment outlet, the differences in discharge, Q, between the free run and each of the 5 "perturbed" simulations were calculated. is the "free run", Q pert is the "perturbed" simulation done with a +10 % perturbation of S (dashed red line) or V (dashed green line); Q is the difference between the "perturbed" simulation and the "free run" for S (in dashed red line with circles) or V (dashed green line with circles).
The drainage coefficient, ds, and the parameter, w, which controls the contribution of drainage to delayed flow, exert a negligible influence on flood discharges. The discharge differences between the "reference" and "perturbed" runs, Q lie below 2 m 3 s −1 , less than 1 % of the peak discharge value in the "free run". The discharge is also relatively insensitive to the K 0 parameter as the maximum difference ( Q) equals 6 m 3 s −1 (approximately 3 % of the free run peak discharge). This parameter has a relatively small impact on the slope of the hydrograph during the rising and recession limbs. As the value of parameter K 0 increases the slopes of the hydrograph flatten. For the transfer velocity, V (Fig. 2), the maximum Q equals 47 m 3 s −1 (approximately 23 % of the free run peak discharge). This parameter influences the flood peak arrival time. As the velocity V increases, the flood peak arrives earlier. The sensitivity of the model to parameter S is also significant (Fig. 2) in that a 10 % perturbation in its value causes a maximum perturbation of 25 m 3 s −1 (approx. 12 % of the "free run" peak discharge). This parameter mainly influences the peak flood intensity. As S grows larger, the catchment's initial moisture deficit is increased, lowering the flood peak intensity. This parameter corresponds to the initial catchment wetness state, with respect to which event-based models are known to be highly sensitive (Zehe and Blöschl, 2004;Berthet et al., 2009). Because simulated discharges are most sensitive to S and V , the data assimilation technique will focus on correcting these two parameters.

Implementation of the assimilation technique
The assimilation algorithm (Sect. 3.1) was implemented for use with the hydrological model (Sect. 2.1) using the Open-PALM dynamic coupling software (Lagarde, 2000;Lagarde et al., 2001), developed at CERFACS (http://www.cerfacs.  fr/globc/PALM WEB/). This software was originally developed for the implementation of data assimilation in oceanography for the MERCATOR project. PALM allows for the coupling of independent code components with a high level of modularity in the data exchanges and treatment while providing a straightforward parallelization environment (Fouilloux and Piacentini, 1999;Buis et al., 2006).
The objective here is to improve flood peak forecasting using a record of past flood events. Observed discharges are assimilated during an assimilation window ("assimilation period") that extends from the beginning of the event until time t 0 , 3 h prior to the flood peak (Fig. 2), in order to correct model parameters. This lead time arbitrarily fixed 3 h before the peak represented a compromise between the usefulness of the lead time for forecasters and the frequency of discharge observations available for the data assimilation technique. The corrected parameters are then used to integrate the hydrological model over the "assimilation period" and beyond, until the end of the event ("forecast period"). During the "forecast period", the simulated discharges are calculated with a known future rainfall (observed rainfall is used to force the model in order to "forecast" a past event). This choice allowed for the assessment of the performance of the data assimilation method without being masked by uncertainties in the forecasted rainfall. However, this methodology did not allow us to test the performance of the hydrological model coupled with the data assimilation method in a realtime framework which is beyond the scope of this paper.
Following from the results of the sensitivity study, the control vector x = (S, V ) T contains the two most influential parameters S and V . Their a priori values are stored in the background vector, x b = (S b , V b ) T . S b is given by the linear regression with Hu2 (Sect. 2.2) and V b is the transfer velocity obtained by the global calibration over the 21 events. The background vector is used to compute the background simulated hydrograph. The standard deviation of the error in S b was fixed at 19 % of S b . This percentage corresponds to the ratio of the average linear regression residuals to the average of the calibrated S values. To estimate this percentage, the mean deviation between Sreg and Scal was estimated and then divided by the mean value of Scal. The background vector is used to compute the background simulated hydrograph. The standard deviation of the error in S b was fixed at 19 % of S b . This percentage corresponds to the ratio of the average linear regression residuals to the average of the calibrated S values. The standard deviation of the error in V b was set to 0.2 m s −1 , which is the standard deviation value obtained following an event-by-event calibration. The errors on S b and V b are assumed to be uncorrelated. In the following experiments, data assimilation is tested with a control vector that contains (i) only S, (ii) only V or (iii) both S and V .
The observations used and stored in the observation vector y o correspond to the first n observed discharges since the beginning of an event at the catchment outlet. The errors in the observed discharges are also assumed to be uncorrelated, leading to a diagonal matrix R. The observation error standard deviation is set to 20 m 3 .s −1 for discharges lying between 20 and 300 m 3 s −1 (Sect. 3.4). Observations above 300 m 3 s −1 are not assimilated due to potentially significant error resulting from the extrapolation of the rating curve. Observations below 20 m 3 s −1 are not assimilated due to significant error in this range since the hydrological model was calibrated using flood discharges.
The observation operator H is the hydrological model forced by rainfall inputs. Applied to x this operator produces a record of simulated discharges y = H (x). The computation of the Jacobian matrix, H l (associated with the non-linear observation operator H ), using a finite difference scheme, requires several runs of the hydrological model, namely: -One run with the reference parameters x l .
-An additional run for each perturbed parameter: x l + dS, 0) T for the perturbation applied to S, and x l + 0, dV ) T for the perturbation applied to V .
The dS and dV values were chosen to be of the same order of magnitude as the difference between the background and the analysis. Furthermore, the hydrological model (i.e. H ) remains almost linear over the dS and dV intervals. To guarantee that the linear assumption is respected, we check that H (S + dS) and its linear approximation H (S) + H dS are nearly the same. For events where the correction (i.e. the interval between the background and the analysis parameter values) is notable, an outer loop is used. Using the Open-PALM coupler, the computation cost of calculating the Jacobian was limited by running these model runs in parallel. When the discharge observation (Fig. 2, blue cross) exceeds the corresponding discharge simulated with the background control vector (Fig. 2, black curve), the correction estimated by the data assimilation algorithm tends to decrease the initial deficit of the soil moisture reservoir in order to increase the analysis discharge simulated over the assimilation  period (Fig. 2, red curve). Since the analysis parameters are used to compute the discharge over both the assimilation and forecast periods, the correction of S leads to a monotonic correction of the discharge over the flood event. Since the background discharges are underestimated over the entire event, the assimilation of one observation allows for an improvement of the entire flood simulation, especially the peak discharge. For the background run, the peak discharge was underestimated by 16 %, whereas after assimilation the peak is overestimated by 7 %. The error in the flood peak estimation was thus reduced by 9 % compared to the observed value. If the difference between the background simulated discharges and the observations is not monotonic, the correction calculated by data assimilation is insufficient (when the sign of the difference is constant) or, worse, leads to a degradation of the discharge simulation (when the sign of the difference changes during the flood event). This last case is discussed in Sect. 4. Figure 3 provides an illustration of how the outer loop operates. For each iteration of the outer loop, the non-quadratic cost function J (x) (black curve) is approximated at the reference point x l (blue cross for the first iteration, red for the second) by a quadratic incremental cost function J inc l (δx) (blue curve for the first iteration, red for the second), which has the same gradient as J (x) at the reference point x l . The initial approximation of J (x) by J inc 1 (δx) is calculated about the background x 1 = x b (blue cross). Next, the BLUE algorithm gives the minimum, x a 1 (blue circle) of J inc 1 (δx). This minimum serves as the linearization point x 2 = x a 1 (red cross) during the next outer-loop iteration. As seen in Fig. 3, x a outer-loop iterations were performed due to computation cost constraints. Nevertheless, this number of iterations still assures that the minimum of the incremental cost function J inc l (δx) (red and blue circles) converges to the cost function minimum J (x) (black circle). In our case the computational cost constraints with 5 iterations are low (few minutes) in comparison with the observation frequency (1 observation/h). In the case of a sliding time window, many more model integrations would be necessary and one would want to reduce the number of iterations to three.

Efficiency of the data assimilation method
To assess the efficiency of the data assimilation method in forecasting the peak flow, the following criterion is defined: where Q obs p is the observed peak discharge, and Q sim p the simulated peak discharge. E Qp denotes the relative deviation with respect to the observed peak discharge; it can be calculated either before assimilation (E b Qp computed with the background peak discharge, Q b p ) or after (E a Q computed with the analysis peak discharge, Q a p ). If E Qp = E a Qp − E b Qp < 0, then the assimilation procedure has improved the peak discharge simulation.
The criterion E tp allows for the evaluation of the impact of the data assimilation procedure on the peak discharge arrival time and is defined by where t obs p is the observed peak discharge arrival time and t sim p the simulated peak discharge arrival time. This offset can be calculated before assimilation (E b tp computed using the background peak discharge arrival time t b p ) or after (E a tp computed using the analysis peak discharge arrival time t a p ). If E tp = E a tp − E b tp < 0, then the assimilation technique has reduced the time offset between the simulated and observed peaks.

Statistical interpretation of the S and/or V correction
The assimilation procedure was applied to 12 flood events containing 20 flood peaks in order to correct: (i) parameter S alone, (ii) the transfer velocity, V alone, and iii) both S and V simultaneously. Table 2 presents the value of the parameters S and V before (background value) and after (analysis value) correction by the data assimilation method.
As expected, the S correction modifies the flood peak intensity but not the peak arrival time as show in Fig. 4a and d. The assimilation improves the flood peak estimation by an average of 12 % over 14 of the simulated peaks (events positive values in Fig. 4a). For 2 peaks (events with no bar in Fig. 4a), there are no observations above 20 m 3 s −1 during the assimilation period, thus no data points are assimilated. For 4 peaks (events with negative values in Fig. 4a), the assimilation procedure leads to the degradation of the peak simulation for reasons that will be discussed in Part 5.
The parameter V correction has a less pronounced effect on the flood peak intensity (Fig. 4b and e). For 9 of the 20 peaks tested, the V correction improved the peak intensity by 8 % on average (negative E Qp in Fig. 4b). For the other 9 peaks, the peak intensity was unchanged following the V correction (zero E Qp in Fig. 4b), and for the last 2 peaks the correction degraded the quality of the flood peak estimation (positive E Qp in Fig. 4b). As expected, the correction of this parameter also modified the flood peak arrival time (non-zero E tp ). For 13 of the 20 peaks tested, the peak arrival time remained unchanged (zero E tp in Fig. 4e). For 2 of the peaks, the time offset between the simulated and observed peaks was reduced by 1 h (negative E tp in Fig. 4e), and for the 5 remaining peaks, the offset increased by 1 h (positive E tp in Fig. 4e). The correction of V by assimilating discharges at the beginning of the flood tended to degrade rather than improve the simulated peak arrival time. This is because flood discharges at the start of the event correspond to the arrival of runoff located near the catchment outlet. The transfer effect is thus limited for these initial discharges, and the difference between simulated and observed discharges stems from the runoff production. Because it is not the source of uncertainty in this case, correcting a transfer function parameter with the initial flood discharges may introduce major errors into the discharge forecast. A correction of the value of the threshold that triggers the direct runoff could be more efficient. This threshold has an influence on the first discharges of the rising limb.
Correcting S and V simultaneously allows for the modification of both the flood peak intensity and arrival time as presented in Fig. 4c and f. For 14 of the 20 peaks, the correction of both S and V serves to improve the flood peak estimation by an average of 14 % (negative E Qp in Fig. 4c). For 2 peaks the correction had no impact (zero E Qp in Fig. 4c). For 4 peaks the S and V corrections deteriorated the estimation of the peak onset (positive E Qp in Fig. 4c). The flood peak arrival time E tp was also affected. In the majority of cases (for 5 events out of the 6 with an altered arrival time following assimilation), the time offset between the simulated and observed peaks increased (positive E tp in Fig. 4f), as presented in Fig. 4e.
To summarize the results presented in Fig. 4, in more than 75 % of cases, regardless of the corrected parameter, data assimilation improves the estimation of the peak discharge indicated by the negative value of E Qp . Nevertheless, the correction of V on its own or accompanied by the S correction actually hinders simulation quality, increasing the offset between the simulated and observed peak arrival times. It is thus preferable to avoid correcting V and to correct S alone.

Sensitivity to the ratio of the matrix B to the matrix R
The efficiency of the assimilation procedure depends on the ratio of the matrix B to the matrix R. The data assimilation procedure was applied while varying the observation error standard deviation between 0.01 and 100 m 3 .s −1 and keeping the background error standard deviation equal to 19 % of the S component of x b . This experiment aims at demonstrating the validity of the data assimilation procedure as well as confirming that the value selected for the observation error standard deviation σ R in this study is in agreement with a physically reasonable value. σ R indicates the relative error in the river discharge observations. Figure 5 shows that, as expected, data assimilation has a larger impact when the observation error standard deviation is smaller, meaning that the analysis is significantly different from the background since more confidence is given to the observations. For a standard deviation of 1 m 3 s −1 and below, the observations are assumed to be almost perfect and  When the observation error standard deviation σ R increases above 1 m 3 s −1 , the confidence given to the observations is reduced, and the effect of the assimilation algorithm gradually decreases. When σ R increases, the boxplot median tends toward zero and the boxplot spread becomes narrower. While the third quartile is equal to zero for σ R greater than 10 m 3 s −1 , the first quartile becomes less and less negative.
When σ R is equal to 10 m 3 s −1 , the flood peak simulation is degraded for the October 2001 and September 2005 peaks, which are always both positive outliers. The deterioration produced by data assimilation can be explained by an underestimation of the observation error standard deviation for October 2001 and by a difference in the rainfall estimation between the 2 peaks for September 2005. For a standard deviation of 20 m 3 s −1 , only the September 2005 peak is considerably degraded, which represents the only positive outlier (red cross). The third quartile value is lower than for 10 m 3 s −1 , while the first quartile value remains high (relative to the other tests). For these reasons this standard deviation value (20 m 3 s −1 ) was chosen for data assimilation.

Discussion
The correction of S leads to a monotonic correction of the flood discharge over the flood event. If the difference between the background simulated discharges and the observations is not monotonic, the correction calculated by data assimilation is negligible (when the sign of the difference is constant) or, in the worst case, leads to a degradation of the discharge simulation (when the sign of the difference changes between the assimilation and forecast period). This section details situations in which the data assimilation procedure leads to a degradation of the flood peak simulation as observed for the following events: the October 1994 (1st peak), December 1997, December 2002, and September 2005 (2nd peak). Two situations were identified as leading to a deterioration in estimation quality: (i) the model does not reproduce the rising limb of the hydrograph as quickly as it occurs in the observations (a situation encountered in all four problem cases); and (ii) in the case of floods with multiple peaks (i.e. the fourth December 2002 peak and the second September 2005 peak), model error in the peak discharge estimation differs from one peak to the next.
The purpose of this discussion section is to analyze these two situations, provide a solution to improve results and to develop a set of limits for the use of this data assimilation approach.

Model error in the simulation of the rising limb
For the four study peaks which were degraded following assimilation, the model did not reproduce the flood rise as quickly as the observations. Consequently, the sign of Figure 6. Data assimilation results for the December 1997 event. The blue curve corresponds to observations; the black curve is the background simulated hydrograph (before assimilation), and the red curve is the simulated analysis hydrograph (after assimilation). Blue crosses represent the assimilated data points. Fig. 6. Data assimilation results for the December 1997 event. The blue curve corresponds to observations; the black curve is the background simulated hydrograph (before assimilation), and the red curve is the simulated analysis hydrograph (after assimilation). Blue crosses represent the assimilated data points. the difference between the background and observed discharges changes during these events. As seen with the December 1997 event (Fig. 6), the model without assimilation (black curve) underestimates the discharge at the beginning of the flood rise (before t inter), while it overestimates the discharge at the end of the rise (after t inter). This situation can be deconstructed by noting that during the flood rising phase, the observed and simulated hydrographs intersect. Since the assimilation period includes a large number of discharges underestimated by the model, the assimilation technique tends to compensate for this underestimation by lowering the S value. This correction exacerbates the peak discharge overestimation (see Fig. 6, red curve: hydrograph after assimilation). In these 4 problem cases, an assimilation that uses observations well before the peak discharge (i.e. prior to the intersection of the rising limbs of the hydrographs) uses excessively high discharge observations in comparison with the simulated discharges that serve to amplify the flood peak overestimation.
In order to improve the assimilation results for these test cases, we increase the value of the threshold above which observed discharges are assimilated so as to only assimilate discharges positioned after the intersection of the rising limbs. Hence, this threshold has been raised from 20 to 60 m 3 s −1 ( Table 3).
The threshold increase from 20 to 60 m 3 s −1 limits the deterioration of the fourth December 2002 peak ( E Qp , decreases from +0.17 for the threshold at 20 m 3 s −1 to +0.13 for the threshold at 60 m 3 s −1 ) and improves the December 1997 peak forecast ( E Qp , decreases from +0.06 for the threshold at 20 m 3 s −1 to −0.04 for 60 m 3 s −1 ). However, as the threshold is increased there are no more observations left to assimilate for the October 1994 (1st peak) event. The case of September 2005 (2nd peak) will be discussed in Sect. 5.2. In order to adjust the observation threshold, the observed hydrograph must be described up to the flood peak, which limits the use of this solution to past events only. Please note that such an increase in the threshold value is acceptable as a sensitivity test but not for real-time applications. If this higher threshold limits the effects of the data assimilation method for these four problematic peaks, the results would be probably less positive for all other peaks. An alternative approach compatible with real-time forecasting would be to carry out the data assimilation analysis using a sliding time window, including observations as they become available to correct the catchment wetness state over time. This sequential approach assumes that the uncertainty results not only from the estimation of the initial condition of the model given by the wetness state indicators, but also stems from the evolution of the wetness state given by the model equations during the event. With such an implementation, each observation would be used once. At each time step, the analysis from the previous window becomes the background for the current window and the B matrix is evolved in time by the Kalman filter equations.

Challenges in representing multiple peak episodes presented by a sample case
Among the four peaks for which assimilation degraded simulation quality, two of them occurred following an initial discharge peak. The loss of model quality observed after data assimilation may thus be explained not only by a problem in representing the flood rise, but also by an error in the estimation of the rainfall peak, which differs from one peak to the next. Such is the case for the second September 2005 peak which displays 2 successive flood peaks (Fig. 7). The model without assimilation (black curve) underestimated the first observed flood peak (blue curve), while overestimating the second. If the error is not monotonic over the 2 peaks, then assimilating the discharge data from the first peak in order to correct the second could lead to significant errors. To identify the consequences of errors resulting from data assimilation, the discharge data were then assimilated in two ways: (i) data from both the first discharge peak and the beginning of the second peak were assimilated ("grouped assimilation"); and (ii) only discharge data from the beginning of the second peak were assimilated ("separate" assimilation). In the first  The blue curve corresponds to observations; the black curve is the background simulated hydrograph (before assimilation), and the red curve is the analysis simulated hydrograph (after assimilation). Blue crosses represent the assimilated data points.
case, the discharge data from all peaks were assimilated together in a "grouped" manner, while in the second case the data from each peak were assimilated "separately". The separate assimilation was tested on both the fourth December 2002 peak and second September 2005 peak. For the latter, "separate" assimilation still deteriorated the flood peak estimation ( E Qp = +0.17), albeit to a lesser extent than the grouped assimilation ( E Qp = +0.42). For the fourth December 2002 peak, the separate assimilation approach improved the flood peak estimation ( E Qp = −0.14), whereas the grouped approach still degraded the flood peak estimation E Qp = +0.17). The fact that the assimilation is unable to improve both peaks reveals that model errors do not result solely from an inaccurate specification of S. Other sources of errors have to be considered and corrected by the data assimilation technique such as the threshold that triggers the runoff or the rainfall that force the model. The correction of this last component is within the scope of Harader et al. (2012). In the case of the second or subsequent peaks in multi-peak events, such as December 2002 or September 2005, uncertainties in the parameters or variables which control the shape of the falling limb, such as ds and w or stoc(t), could be also taken into account.

Conclusions
Hydrological model calibration is essential for ensuring that simulations are coherent with observations. The global calibration approach is insufficient as certain model parameters (as well as their errors) are event dependent. Within the framework of data assimilation, these parameters can be calibrated using observational data as they become available in order to forecast the flood peak. In this study flood forecasting is done using known future rainfall, while acknowledging that this would not be the case in a real-time framework.
The assimilation technique presented in this paper uses discharges observed at the catchment outlet during the beginning of the flood event in order to correct the initial deficit of the soil moisture reservoir, S and/or the transfer speed V , which were found to be the most sensitive model parameters. As a consequence, in most cases, the use of the corrected parameters allows for an improved representation and forecast of the entire flood event.
In this study a BLUE algorithm together with an outer loop procedure were built around a distributed, event-based, parsimonious hydrological model. The outer loop accounts for some of the non-linearities between the model parameters and simulated discharges. The assimilation experiments conducted in this study demonstrated that the correction of S is preferred over the correction of V . The S parameter influences the production function which plays an essential role in the representation of discharges early in the flood event (during the assimilation period), leading to an improvement in the flood peak estimation for 75 % of the tested peaks when correcting S. In contrast, the error in the observed discharge at the beginning of the event is not related to errors in the transfer function since its effect on the simulated hydrograph is limited during this period; thus, the correction on V is not recommended for this data assimilation approach.
Given that model parameters are assumed to be constant over the flood event, the correction of the discharges introduced by the algorithm is monotonic. Such a correction is not well adapted to cases in which the difference between simulated and observed discharges strongly varies over the flood event, as illustrated for 25 % of the test cases. These results can be improved by either using only observations above an arbitrary threshold or separating multiple peak events into distinct events. Possible improvement could be expected by extending the data assimilation correction to other components of the hydrological model such as the value of the threshold that triggers the direct runoff, the parameters that control the shape of the falling limb during multi-peak events or the rainfall calculated over the catchment. This last correction is addressed in Harader et al. (2012). The present study, which presents the benefits and the limitations of a data assimilation for hydrological forecast, must be extended to other lead times, other catchments and hydrological models in order to draw more general conclusions. A potential follow-up to this study is the sequential application of the data assimilation algorithm on a sliding time window as in Ricci et al. (2011). This technique would take into account the latest observations available to provide a time varying estimation of the model parameters over the flood event. This sequential approach would deal with two kinds of uncertainties: the estimation of the initial condition of the model given by the indicators at the beginning of the event, and the evolution of the catchment wetness state given by the model during the event. This kind of sequential assimilation could be adapted to a real-time application.