Interactive comment on “ Correcting the radar rainfall forcing of a hydrological model with data assimilation : application to flood forecasting in the Lez Catchment in Southern France ”

General comments The authors use a distribute rainfall-runoff model to reconstruct flood events in a watershed of southern France (Lez Catchment). Their work focus on the correction of radar data assimilation. The paper is good in its contents, this is the reason why I would recommend publication, but, on the other side, it is in many of its parts wordy, whereas many other parts should be extended. The improvements on this side can for sure put more in evidence the good scientific content of the paper. I hope that my comments will be useful.


Introduction
For flash flood prediction, hydrologists may use tools as rudimentary as rainfall-discharge curves or as refined as complicated physical and distributed hydrological models, all with the goal of converting atmospheric and soil conditions into discharge volumes, flood peak amplitudes and arrival times.All of these tools are subject to uncertainties related to their inputs and parameterisations.Rainfall-runoff models are sensitive to rainfall quantities and their spatial distribution throughout the catchment, as runoff generation depends upon rainfall location.Errors in rainfall estimates have a significant impact on prevision and event reconstruction quality.

E. Harader et al.: Data assimilation for the correction of radar rainfall
In studies of flash flood modelling for Romanian catchments between 36 and 167 km 2 , Zoccatelli et al. (2010) demonstrated that neglecting the spatial variability of rainfall resulted in a deterioration of the simulation quality.Roux et al. (2011) showed that the MARINE model (Estupina-Borrell, 2004) was dependent upon the distribution of rainfall data in order to correctly represent the soil saturation dynamics of the 545 km 2 Gardon d'Anduze catchment in Southern France.The sensitivity of models to rainfall distribution highlights the importance of using a rainfall product with a fine spatial and temporal resolution, such as that provided by weather radar.
However, the use of radar data is often limited by increased uncertainties compared to ground rainfall measurements due to nonlinearities in the rainfall-reflectivity relationship, ground clutter and beam blocking (Borga, 2002).In the Cévennes region where the Lez catchment is located, a hilly terrain complicates the process of separating rainfall and terrain backscatter.Pellarin et al. (2002) demonstrated that selecting the scan used in mountainous regions based on distance considerations, as done so for the HYDRAM rainfall product used in this study, leads to a lower quality rainfall product compared to using a composite (highest quality scan at any given point) method (Cheze and Helloco, 1999).Additionally, in the Lez catchment, radar data quality varies by season and is diminished in winter months due to bright band effects related to predominantly stratiform rainfall (Coustau et al., 2011;Emmanuel et al., 2012;Tabary, 2007).A possible post-treatment correction to radar rainfall is the removal of the mean field bias (MFB) (Wilson and Brandes, 1979), a correction which uses rain gauge data to eliminate errors due to instrumentation and a nonlinear vertical profile reflectivity (VPR).Adjustment of radar rainfall using rain gauge data has been shown to lead to improved prediction accuracy (Vieux and Bedient, 2004;Cole and Moore, 2008).
Identifying a correction to the rainfall data input to hydrological models can also be formulated as an inverse problem (Tarantola, 2005;McLaughlin and Townley, 1996) solved in the framework of data assimilation.Data assimilation (DA) for the improvement of hydrological event reconstruction or forecast has been already demonstrated as effective, e.g., Aubert et al. (2003), Moradkhani et al. (2005), Pauwels et al. (2001), Thirel et al. (2010) and Vrugt et al. (2005).However, previous literature has often focused on correcting rainfall without the direct implication of downstream hydrological applications.Chumchean et al. (2006) used a Kalman filtering approach, modelling the logarithmic MFB as an autoregressive process.Seo et al. (1999) recursively calculated the MFB using an exponential smoother.These techniques are convenient in that they depend only on radar and rain gauge measurements, however, studies of their impact on hydrological modelling could be further developed.In the context of flood forecasting, Kahl and Nachtnebel (2008) adopted an updating technique which relates the rainfall correction to a hydrological simulation through the minimisation of an objective function.However, the objective function has two drawbacks: (i) it has no explicit solution and (ii) it does not take into account the observation error.This study builds upon established methods by using DA to correct rainfall while focusing on downstream hydrological applications.
A common approach to data assimilation is the Kalman Filter (KF) algorithm.The KF corrects a set of a priori parameters and/or model states (the background) stored in the control vector using observations to produce a set of optimal model states or parameters (the analysis).Assuming that the observation operator mapping the control vector onto the observation space is linear, the algorithm calculates the analysis by performing a linear combination of the background and analysis, each weighted by their respective error covariances.The extension of the KF to nonlinear operators (Extended Kalman Filter -EKF) implies the computation of a local estimate of the tangent linear of the observation operator (Goegebeur and Pauwels, 2007).The EKF analysis is similar to the incremental 4D-Var (4D-inc) analysis in that they both rely on the local linearisation of the observation operator (Bouttier and Courtier, 2002).Both variational and filtering analyses are based on the minimisation of a cost function that describes discrepancies between simulated and observed values as well as their associated error statistics.However, these algorithms differ in the way the minimisation is performed: variational techniques use a minimiser and are adapted to large dimension problems, whereas filtering techniques explicitly solve for the analysis using matrix multiplication and inversion that are only affordable for small dimension problems such as the one presented here.The limitations of both 4D-inc and EKF are due to the use of a local estimation of the tangent linear of the observation operator and can be partly overcome with an update of the linearised operator also called an outer loop (Thirel et al., 2010).Another possible alternative to the EKF, the Ensemble Kalman Filter (EnKF) estimates error statistics from an ensemble of model runs and enables a stochastic estimate of the covariance matrices taking into account the nonlinearity of the observation operator (Weerts and El Serafy, 2006;Pauwels and De Lannoy, 2009;Moradkhani et al., 2005).Ensemble methods, such as the EnKF, Particle Filter or the Maximum Likelihood Ensemble Filter (MLEF), can thus be used for nonlinear systems; however, the quality of the resulting analysis strongly depends on the initial sample and whether it does or does not properly represent the uncertainty of the system.
The main objectives of this study are: (i) to assimilate discharge data using an EKF to correct radar rainfall data which is a key source of uncertainty in hydrological modelling and (ii) to apply this correction to flood simulation and forecast in order to examine the quality of hydrological prediction using DA.Other uncertainty also exists in the model structure and physical catchment properties, but radar rainfall was selected as the target of DA because it is a key factor in the hydrology of the catchment and it provides several advantages over rain gauge data if its uncertainty can be reduced.In order to evaluate the quality of the DA correction, comparisons were made with the MFB; then, to evaluate the predictive capacity of the method, the correction was applied in a forecast-like setting.
The paper is outlined as follows: Sect. 2 includes a description of the study site, the model structure and calibration, the DA procedure, a description of the experimental setup and examples of assimilation performed in reanalysis and pseudo-forecast mode.The results of the study and the impacts of data assimilation on the efficiency of the hydrological model are then presented in Sect.3. Finally, a summary of the key results obtained and conclusions are discussed in Sect. 4.

Study site: the Lez catchment
The Lez catchment in Southern France (Fig. 1) is a mediumsized karstic basin located in the Hérault department, 15 km north of the town of Montpellier.The catchment is 114 km 2 at Lavalette, where discharge measurements are taken.This portion of the Lez river is fed by several upstream tributaries: the Lirou, Yorgues and Terrieu.The Lez River stretches for 26 km between its source and the Mediterranean Sea.
The landscape of the Lez catchment at Lavalette is defined by plains and hilly garrigue with limestone outcrops and very little urbanisation.The plains are composed of 200 to 800 m thick Valanginian marls (a mixture of calcium carbonate and clay minerals formed during the Early Cretaceous period), covered by soil usually less than 1 m thick.Land use ranges between agricultural (vineyards) and forest in the plains, along with undeveloped garrigue; the limestone outcrops have very little soil cover and thin vegetation.
The source of the Lez is a seasonal spring which serves as the main outlet of a 380 km 2 limestone and dolomite karstic aquifer (shown by the dotted line in Fig. 1) (Avias, 1992).Karstic systems are defined by the presence of conduits and fractures in the underlying limestone bedrock, resulting in complex transport networks and variable response times following rainfall events.The subsurface processes that contribute to runoff are poorly understood: they may reduce flood intensity by storing water in the epikarst and through deep infiltration (Dörfliger et al., 2008) or they may intensify the flood severity through the contribution of groundwater to peak flow (Kong A Siou et al., 2011).

Climate and rainfall data
The climate of the region is generally dry, with mean annual potential evapotranspiration (1322 mm at Mauguio for the period from 1996 to 2005 -Fig. 1) greater than mean annual rainfall (909 mm at Prades for the period from 1992 to 2008).Mean annual evapotranspiration was calculated using the Penman-Monteith equation; this calculation is not available at the Prades rain gauge.Most of the yearly rainfall is received in fall and winter in the form of heavy climatic and orographic precipitations.To the North of the Lez catchment, frontal systems are strengthened by relief changes in the Massif Central.Extreme rainfall events, particularly in late summer and fall periods, are favoured in this region due to humidity generated by the warm Mediterranean Sea and a closed cyclone which helps to transport warm, moist air masses to the coast (Nuissier et al., 2008).In September of 2002, rainfall totalled as much as 600-700 mm over a 24 h period in certain regions (Boudevillain et al., 2011).
Rainfall in the Lez catchment is measured by both an Sband radar located in Nîmes at a distance of approximately 65 km from the basin and a network of 4 rain gauges (Prades, Montpellier-ENSAM, Maugio, Saint Martin de Londres -Fig.1).Radar data were treated using the HYDRAM algorithm developed by Météo-France (Cheze and Helloco, 1999) for the correction of ground clutters, the vertical profile of reflectivity and the conversion of reflectivity to rainfall using the Marshall-Palmer relationship, where Z is the reflectivity in mm 6 m −3 , R is the radar rainfall intensity in mm h −1 and 200 and 1.6 are empirical constants derived from the drop size distribution.For the HYDRAM treatment, the same Z-R relationship is used for stratiform and convective rainfall (Tabary, 2007).The Nîmes radar produces scans at three different elevations at 5 min intervals: 2.5 • (0-22 km), 1.3 • (22-80 km) and 0.6 • (distances beyond 80 km).These three scans are used to produce a radar image which describes rainfall for areas at different distances to the radar.The lowest unobstructed scan is selected for a given distance range.For the Lez catchment, the 1.3 • scan was used (Bouilloud et al., 2010) to produce cumulative rainfall depths at a spatial resolution of 1 km 2 and a time step of 5 min.A network of 20 rain gauges within a 50 km range of the catchment provided cumulative rainfall data for adjustments using the MFB (Fig. 1), a measure of the ratio of radar to rain gauge rainfall during a specified time period (here the length of the flood event): where G i is the rain gauge measurement at location i in mm, R i is the radar measurement at the same location in mm and n is the number of rain gauges selected.The value of the radar measurement at the gauge location was selected to be the average of the central pixel and its 8 nearest neighbours.The ratio of rain gauge to radar measurements is expected to be greater than 1 for distances between 15 and 80 km from the radar where masking effects play an important role (Cheze and Helloco, 1999).

Rainfall events
Table 1 displays the 19 rainfall events measured by HYDRAM-treated radar for the Lez catchment together with their associated MFB values and peak discharges.In general, events lasted several days and cumulative rainfall was sampled at a time step of 1 h.The episode MFBs were between 0.87 and 1.80, indicating that radar was never more than 45 % away from the "true" rainfall value (assuming absolute confidence in ground measurements) with the exception of November 1997.The very high MFB for this event indicates that either the rain gauges, the radar or both may have not been functioning properly.With the exception of November 2008, all events have MFB values greater than 1, with an average of 1.39.As mentioned in Sect.2.1.1,these values are a feature of the distance between the Nîmes radar and the watershed.Rainfall events were separated into two classes based on their peak discharges: regular events which have a peak discharge greater than 40 m 3 s −1 and very small events which have a peak discharge less than or equal to 40 m 3 s −1 .This classification is used to determine the range of discharges that will be assimilated as discussed in Sect.2.4.

The hydrological model
The hydrological model is event-based, parsimonious and distributed.It operates on independent grid cells with an hourly time step using a derived SCS runoff production function (Gaume et al., 2004) and a Lag and Route transfer function (Tramblay et al., 2011).The calibration and adaption of this model to the Lez catchment are presented in Coustau et al. (2012).

The runoff production function
The runoff production function is the link between the precipitation falling over the catchment and the discharge emitted to surface waters.Not all rain becomes discharge and processes such as infiltration, evapotranspiration, percolation and interception determine the eventual fate of incident rainfall.The SCS method for predicting runoff has been validated for medium-sized watersheds in recent literature (Abon et al., 2011;Han et al., 2012).The ATHYS software, developed by HydroSciences Montpellier (www.athys-soft.org),was used to implement a derived version of the SCS equations for this study (Gaume et al., 2004), where i e (t) is the instantaneous runoff rate (or runoff intensity) with units of mm s −1 , i b (t) is the rainfall rate (or rainfall intensity) in mm s −1 and C(t) is the fraction of rainfall contributing to runoff.C(t) is defined as follows, where P b is the cumulative rainfall depth at time t in mm and S is the potential storage depth of the watershed at the start of the event (potential maximum retention) in mm.
To represent the ability of the soil to regain part of its absorption potential during pauses in the rainfall, this version of the SCS method allows the soil to drain.The volume of water lost to drainage is a function of two conceptual reservoirs: the cumulative rainfall reservoir, level P b (t), and the soil reservoir, level stoc(t), in mm, shown in Fig. 2. The cumulative rainfall reservoir represents the total rainfall depth received and is used to calculate of the portion of the incident rainfall contributing to runoff during the event.The soil reservoir represents the amount of rainfall stored in the soil.A portion of the water lost by this reservoir becomes delayed runoff.The rate of drainage of the cumulative rainfall reservoir and the soil reservoir is described by:  (Bouvier and Delclaux, 1996).
where ds is the drainage coefficient in d −1 .This coefficient represents the removal of water through deep infiltration and evapotranspiration during the event.The drainage coefficients of the cumulative rainfall reservoir and the soil reservoir were selected to be the same.The water lost to the system by the drainage coefficient is considered to be either lost to deep infiltration or to re-emerge as delayed surface runoff, i d (t), calculated by where w is the critical soil depth in mm and S is the same as that appearing in Eq. ( 4).The ratio between S and the critical soil depth determines the fraction of drainage that becomes delayed runoff.As S approaches w (going from high S to low S), the proportion of runoff lost to deep infiltration is diminished and a greater portion of the soil reservoir drainage becomes available as delayed discharge.The critical soil depth was added by Coustau et al. (2012) in order to adapt the SCS equations to the behaviour of karstic watersheds and to ensure the proper behaviour of the watershed during the descending limb of the hydrograph by including the participation of subsurface flows.The delayed surface runoff is then added back to the instantaneous runoff rate to produce the total runoff, i t (t).

The transfer function
Supposing that the production function has created runoff at a certain grid location, this runoff must then be transferred ) is based on a unit hydrograph approach in which the discharge produced by each cell is assumed to follow the form of a decaying exponential.In this way, it is similar to the impulse solution of the kinematic wave approach.However, in the present case, the form of the hydrograph is assumed and imposed upon the runoff generated by each cell.This runoff is independent and does not interact with that of the other cells (Olivera and Maidment, 1999;Maidment et al., 1996).Independent grid cells may be a strong simplification; however, runoff is rapidly concentrated, leading to little or no infiltration or storage during flow routing.This is in contrast the kinematic wave approximation (a simplification of the Saint-Venant equations for shallow water flow) or the Manning equations for open channel flow (Bates and De Roo, 2000).In these cases, the discharges from different cells are allowed to interact and the flow rate will depend upon the depth of the runoff contained within the cell.Despite its simplicity, the Lag and Route function has been shown to perform as well as the Saint-Venant equations for certain cases (Lhomme et al., 2004).The use of independent grid cells with a Lag and Route transfer function was selected because it does not require prior knowledge of the hydrodynamic features of the catchment such as roughness coefficients or hydraulic conductivity and has relatively few parameters to calibrate.
The two parameters which describe the Lag and Route function are: V 0 , the speed of propagation in m s −1 , and K 0 , a dimensionless coefficient used to calculate the diffusion time.The propagation time to the outlet, T m , in seconds describes the lag between runoff production at time t 0 and the arrival of an associated elementary hydrograph at the watershed outlet.It is equal to l m , the length of the flow path from the cell to the outlet in m -calculated using a method of steepest descent in order to produce drainage paths for each cell, divided by V 0 .From the propagation time, the diffusion time K m in s is calculated as the product of K 0 and T m .This coefficient represents the velocity distribution of the runoff as it is transferred from the cell to the outlet.For each grid cell, the diffusion time and propagation time are then used to produce an elementary hydrograph, q(t) in m 3 s −1 , produced by the total runoff i t (t 0 ): where A is the area of the grid cell in m 2 .
To measure the quality of the simulations performed by the hydrological model, the Nash-Sutcliffe efficiency criterion (NS) was selected (Nash and Sutcliffe, 1970).This criterion can be expressed as a function of the error between the model discharge at time j (Q sim,j in m 3 s −1 ) and measured discharge at time j (Q obs,j in m 3 s −1 ), summed over j , squared and normalised by the variance of the measured discharge (σ 2 obs ): where j varies from 1 to N, the total number of observations available for the event.For this study, NS is calculated over the entire length of the rainfall event, regardless of the number of observations assimilated.The window of observations selected for assimilation will be discussed in detail in Sect.2.3.
A second measure of quality is the normalised difference in peak flow between the simulation (Q sim, peak ) and the observations (Q obs,peak ), PH:

Sensitivity of the model to rainfall inputs
In this section, the choice of radar rainfall as the target of data assimilation will be explained and the relationship between discharge and rainfall explored.
Rainfall plays a key role in the estimation of discharges using hydrological models.The model used in this study is sensitive to the quantity and intensity of rainfall and this sensitivity varies depending on previous conditions.As the soil reservoir becomes saturated, a greater proportion of incident rainfall runs off and is emitted as discharge.In this way, the response of the watershed to a linear increase in rainfall is expected to be nonlinear because the behaviour of the soil moisture reservoir after 20 mm of rainfall is not the same after 40 mm of rainfall.To illustrate this phenomena, a linear multiplier of the rainfall intensity, denoted α, was introduced into the model: where i b is the observed radar rainfall rate and i b is the rainfall rate used by the model.Figure 3 displays the discharge as a function of α at 3 h before the flood peak.This time step was selected because it demonstrates saturated behaviour for larger values of α and non-saturated behaviour for small α.
The discharge is highly sensitive to rainfall inputs with values near 0 m 3 s −1 for α = 0.5 and 1000 m 3 s −1 for α = 3.As expected, the relationship is nonlinear.This is due to (i) a nonlinear runoff production function which depends on soil saturation and (ii) the differential equations describing soil and rainfall reservoir drainage.Despite nonlinearities, α was chosen as the target of the DA procedure because of the strong influence of the rainfall input upon model results.

Initialisation and calibration of the model
The hydrological model contains several types of parameters: batch-calibrated parameters, mathematical properties of  the equations and the initial condition of the watershed, S, the potential storage depth of the soil reservoir, which must be calibrated separately for each event.It should be noted that while the language "initialisation" is used here, S is a parameter in this data assimilation system and not a model state, thus, it does not evolve during the event.During the calibration process, a mixture of ground rainfall events and high quality (early autumn) radar rainfall events from 1994-2008 was used in order to minimise the error associated with the parameterisation.The first step was to calibrate ds, a mathematical property of the model equal to the coefficient of the exponential recession limb of the hydrograph.When the rainfall rate is zero, discharge consists entirely of delayed runoff and stoc(t) becomes a decaying exponential with a coefficient of ds.The slope of the semi-log plot of the discharge is then equal to ds, the coefficient of the decaying exponent.This value was determined to be 0.28 d −1 for all events.Next, the batch-calibrated parameters V 0 and w were calibrated by selecting the value which maximises the NS of the simulated discharge for a given event and then averaging over all events.To avoid problems of equifinality (Beven and Freer, 2001) during this step of the calibration process, K 0 was set as a fixed value before calculating V 0 and w.Since the diffusion time K m is a function of both V 0 and K 0 , many values of these two parameters can result in the same velocity distribution at the watershed outlet.The parameters V 0 , w and K 0 were determined to be 1.3 m s −1 , 101 mm and 0.3 (dimensionless), respectively, for all events.
Finally, the initial soil moisture deficit or potential storage depth, represented by the parameter S, must be calculated at the beginning of each event.In reanalysis mode, a posteriori S values, denoted S cal , were calibrated for each rainfall event by maximising the NS of discharge simulations forced with the MFB corrected radar rainfall in order to minimise errors in the parameterisation.In pseudo-forecast mode, the event hydrograph is not known.As a consequence, S must be estimated at the start of the event using known indicators of the catchment wetness state at this time.For example, piezometric readings could be used to estimate the state of the karstic aquifer in the morning if heavy rain was predicted for the evening.In this study, a calibration curve relating S to indicators of the catchment wetness state is used to estimate a priori S values for each episode from measurements of aquifer piezometry or soil moisture indicators derived from surface models (Coustau et al., 2012).These estimated S values are referred to as S reg .
Using the historical record of discharge and rainfall from 1997-2008, calibration curves for S were developed using 3 catchment wetness state indicators: Hu2 (%), the piezometer located at Bois Saint Mathieu (m) and the piezometer located at Claret (m) (Fig. 1).These two piezometers were selected for the quality of their relation to the hydric state of the watershed.The Hu2 indicator is modelled by Météo-France (Quintana-Seguí et al., 2008) and estimates the % soil saturation at the root horizon.The measurements for each event are taken as the value of the indicator at 06:00 a.m.UTC the day of the event.Hu2 data are available for 18 of the 19 rainfall events and piezometer data are available for 14 of the 19 events.
For each indicator, a regression of slope M and y-intercept b was formed using the catchment wetness state indicator as the independent variable and S cal as the dependant variable as shown in Table 2, where S is the parameter described in Eq. ( 4) calibrated for each episode.R 2 is the coefficient of determination for the linear regression between S cal and the physical indicators.To validate each regression, split sample tests were performed.Each regression was performed using only the first half of the data available to construct a "historical period"; the S reg values calculated using the validation regression were then compared with the S reg values calculated using the regression for the entire record.The average and standard deviation of the % difference between these two

E. Harader et al.: Data assimilation for the correction of radar rainfall
S reg values are presented in Table 2.The piezometer at Claret was the most robust indicator during this phase of the validation.A second validation was performed by comparing S reg to S cal during the validation period.The average % difference between S reg and S cal was 0.22, 0.21 and 0.16 for Bois Saint Mathieu, Claret and Hu2, respectively.For this test, Hu2 was the most robust indicator.
The S reg values calculated using the different regressions are shown in Table 3.An analysis of the impact of errors in the parameterisation will be presented in Sect.3.2.1.

Data assimilation methods
A non-sequential EKF with an outer loop was selected for this study.Data assimilation was carried out over a time window which includes several discharge observations assimilated in a single analysis to correct the input rainfall described by weather radar.The control vector is a scalar containing a multiplier of the input rainfall assumed to be constant over a time window which contains the entire flood event.The observation operator mapping the control vector on to the observation space (discharges at the catchment outlet) is represented by the integration of the hydrological model.The linearised version of the hydrological model is calculated locally about a reference value of the control vector using a finite difference scheme.This reference value is initially selected as the background control vector.However, this method is limited by the assumption that the observation operator is linear in the vicinity of the background.To account for nonlinearities in the observation operator, an outer loop was applied to the EKF.The outer loop updates the observation operator using the analysis as the reference value and then calculates a new analysis starting from the background control vector.The main advantages of this algorithm are low computational costs for a small control vector and the simplicity of implementation.Using the EKF described above, DA was carried out for heavy rainfall events occurring within the Lez catchment between 1997 and 2008.The analysis was applied in two modes: reanalysis and pseudoforecast.In reanalysis mode, all available discharge observations during the rainfall event were assimilated.In forecast mode, observations up to 3 h before the peak flow arrival were assimilated in order to reproduce an operational forecasting environment.The resulting rainfall multiplier was then applied until the end of the rainfall episode.This choice of assimilation window is intended to demonstrate the possible performance of the algorithm in a real-time forecasting environment, while acknowledging that the peak arrival time would not be known in this case.
In this application of the EKF, information from the background discharge simulation Q sim,b is combined with observed discharges Q obs to calculate a constant multiplier of radar rainfall inputs, α, which is then used to integrate the hydrological model, producing a corrected discharge simulation as shown in Fig. 4. The rainfall multiplier is calculated over a single time window covering the entire flood event (or until 3 h before the peak flow for the pseudo-forecast mode).As a consequence, this multiplier represents the mean behaviour of the rainfall over each event, as it is constant in time and uniform in space.Discharges simulated by the conceptual hydrological model used in this study have a nonlinear dependence on rainfall inputs.In data assimilation, this relationship can be represented as a nonlinear observation operator H.This operator translates rainfall input i b into discharge data Q sim , using model parameters (such as S and V 0 ) to solve ordinary differential equations for state variables stoc and P b : where x is the control vector containing a multiplicative coefficient of the rainfall intensity, denoted α, presented in Sect.2.2.3 and y is the control vector in the observation space (i.e.discharges).It should be noted that subscripts indicating the time dimension of x and y are not included.This is because x is constant over each rainfall event as previously stated and y gathers together model outputs for each observation time over the rainfall event.The translation of rainfall input to discharges at the catchment outlet is represented in Step 1 of the algorithm schematic diagram (Fig. 5).
Assuming that the errors in the rainfall input and the observations follow a Gaussian distribution, the optimal value of the control vector is the analysis, x a , which minimises the cost or misfit function J (Bouttier and Courtier, 2002):  to produce the analysis model trajectory (Q a in green) during steps 1 through 4. In step 5, the observation operation is re-linearised in the vicinity of the analysis and steps 2 and 3 are repeated to form an "outer loop".
(13) The cost function J is the sum of two terms: (i) the difference between the control vector x and its background value x b and (ii) the difference between the control vector in the observation space and the observation vector y o , weighted respectively by the background and observation error covariance matrices, B and R. The background control vector is selected as x b = (1) (no change to the input rainfall) and the observation vector contains the observed discharges during the assimilation window.
The cost function above is at a minimum when its gradient is null, leading to the expression ∇J (x a ) = 0. To express the cost function gradient, the derivative of the nonlinear observation operator with respect to the control vector is necessary.The Jacobian matrix H of the observation operator H is determined using the Taylor expansion computed around a reference vector x ref , initially chosen as x b (Step 2 in Fig. 5):

E. Harader et al.: Data assimilation for the correction of radar rainfall
Using Eq. ( 14) to nullify the gradient of Eq. ( 13), x a can be determined (Step 3 in Fig. 5): where x a is the EKF analysis, d is the innovation vector, y o − H(x b ) represents the difference between the simulated discharge Q sim,b (Step 1 in Fig. 5) and the observed discharge Q obs , K is the gain matrix, B H T (H B H T + R) −1 , and δx is the increment applied to the background.The hydrological model can now be integrated using the analysis rainfall multiplier stored in x a to provide a new estimate of the simulated discharge Q sim,a (Step 4 in Fig. 5).
The use of the Extended Kalman filter analysis equations relies on the hypothesis that H(x) can be approximated as locally linear in the vicinity of x ref and that this approximation is valid on [x a , x ref ].The innovation added to the background is assumed to be sufficiently small that the residual between H (x b + δx) and H(x b ) + H x b δx is negligible for an increment δx applied to the background.Limitations of the non-sequential EKF occur when the innovation extends outside the region where the linearity assumption holds.To compensate for nonlinearities in H(x), the outer loop procedure (Thirel et al., 2010) in Fig. 5 allows for the recalculation of the linear tangent H at the location of the analysis of the previous iteration x a (Step 5 in Fig. 5) in order to create a new quadratic approximation of J , as shown in Fig. 6.At the optimal value of the analysis, the minimum of the quadratic approximation and the non-quadratic cost function will coincide.By re-calculating the linear tangent about the analysis, the minimum of the quadratic cost function approaches that of the non-quadratic cost function.The analysis calculated with the new quadratic approximation then provides an improved estimate of the non-quadratic cost function minimum.This method could also be applied to a 4D-Var incremental algorithm.
The B matrix represents the background error covariance, which is the error in the rainfall multiplier.This error is assumed to follow a Gaussian model and is described by its variance as the control vector is a scalar.However, the variance of the rainfall corrective coefficient is difficult to define because it is the uncertainty in a correction applied to the radar rainfall and not the uncertainty of the measure itself.In order to define B, α was assumed to have an error near that of the MFB, which has a standard deviation of 30 % and an average deviation of 40 %.The standard deviation of α was selected as the higher of these two error estimates as a precaution.
The observation errors are supposed uncorrelated, making R a diagonal matrix.A proportionality coefficient, β obs was used to calculate the observation error variance σ 2 obs in order to control the amount of confidence placed in observations depending on the assimilation window selected (reanalysis or pseudo-forecast mode): σ 2 obs,i = max where t i is the initial time step and t f is the final time step of the assimilation period.R has a lower bound of 0.01 m 6 s −2 and no upper bound.As the errors coming from each source of information are not precisely known in pseudo-forecast mode, different values of the proportionality coefficient were considered as described in Sect.2.4.In reanalysis mode, β obs is selected such that all discharges above 2.5 m 3 s −1 have the minimum error covariance of 0.01 m 6 s −2 .This choice is based on the use of the Nash-Sutcliffe criterion for measuring model efficiency.The Nash-Sutcliffe criterion measures model outputs against observed data, placing absolute confidence in the observations.In order to improve this criterion, the algorithm was used to match observations as closely as possible in the reanalysis mode.In pseudo-forecast mode and for discharges below 2.5 m 3 s −1 in reanalysis mode, the observation error becomes heteroscedastic (variance changing).
The variance is then proportional to the discharge measurement as in Moradkhani et al. (2005), though an inverse proportionality scheme was selected in this case in order to place more weight on high flow conditions.This is better suited for peak flow forecasting.

Experimental set-up
For the 19 radar rainfall events, the range of assimilated discharges is 15-300 m 3 s −1 for normal episodes and 2-40 m 3 s −1 for very small episodes (peak discharge less than or equal to 40 m 3 s −1 ).Very large discharges are unreliable due to the use of a rating curve to calculate the river stagedischarge relationship beyond 300 m 3 s −1 .Small discharges are eliminated in order to better represent the flood behaviour of the watershed.For each calculation of the analysis control vector in both reanalysis and pseudo-forecast modes, five iterations of the outer loop method were used.Data assimilation was applied to all episodes in both pseudo-forecast and reanalysis mode, with the exclusion of October 2008.The rising limb of this event takes place over a period of time less than three hours long, thus, no discharge measurements are assimilated in pseudo-forecast mode.
Episodes with notable double peaks (September 2002, October 2002, December 2002, September 2005and October 2008) are separated into single peaks prior to assimilation due to the inability of the hydrological model to properly simulate multiple peaks in succession.The model has difficulties in representing the initial state of the catchment at the start of the second flood peak.This may be due to the influence of the karst in sustaining the discharge during the recession limb of the hydrograph (Coustau et al., 2012) or the effect of random variations in the rainfall error.By separating the peaks, data assimilation may help to correct some of the temporal variations in the rainfall error.

Reanalysis mode
In reanalysis mode, the initial deficit of the soil moisture reservoir is parameterised by S cal and β obs is chosen to be 0.25 m 6 s −2 in order to reflect an almost complete confidence in the observations.Results of the reanalysis mode are first compared to the background simulation in Sect.3.1.1and then to simulations forced with MFB-corrected rainfall in Sect.3.1.2.DA was not applied to the simulations used for comparison.
To illustrate the DA procedure, the episode of November 2008 was selected.In reanalysis mode, the potential storage depth of the catchment, S cal is 142 mm.β obs is chosen to be 0.25 m 6 s −2 .As shown in Fig. 7, the NS is improved from −0.52 to 0.72 following assimilation.In this case, α = 0.70 for the analysis, meaning that the optimal state of the rainfall is less than that predicted by the uncorrected radar data.The reduction in the amount of rainfall then results in an analysis hydrograph that is smaller than the background hydrograph.

Pseudo-forecast mode
Several modifications to the assimilation procedure are necessary to assimilate data in pseudo-forecast mode.First, the observation error covariance, parameterised by β obs , must be adjusted to reflect representativeness errors due to a reduced number of observations being assimilated (only the start of the event is known).It is expected that β obs will need to be increased in this case to reflect less confidence being placed in the observations.Next, an a priori estimation of S (S reg ), as presented in Sect.2.2.4, is required.The lack of a fully-described hydrograph leads to uncertainties in pseudo-forecast mode that are not present in the reanalysis.In order to characterise the uncertainty in the observations, an initial experiment was carried out by assimilating discharge data using S cal and three different values of β obs .S cal was used so that parameterisation errors would not influence the results.These results were then compared to assimilation using S hu2 as presented in Sect.3.2.1.Using the β obs determined in Sect.3.2.1,experiments using the 3 different S reg parameterisations, Hu2, Bois Saint Mathieu and Claret, are presented in Sect.3.2.2.The goal of this test is to characterise the impact of the parameterisation upon the results and to determine if certain catchment moisture state indicators provide better S values than others.The experiments are measured against assimilation using S cal , which should have the best performance due to an improved parameterisation.Simulations using the Hu2 parameterisation are expected to have the lowest performance, since this catchment wetness state indicator contains model error in addition to measurement uncertainty.
In pseudo-forecast mode, observations are assimilated from the start of the event until 3 h before the peak discharge.This process is illustrated in Fig. 8 for November 2008; the first and final analyses of the outer loop are shown.For this demonstration, S and β obs were kept the same as those for the reanalysis mode.The first iterate of the outer loop has the best NS with 0.71 which is nearly equal to that of the reanalysis mode.However, this is not the optimal state for the assimilation period (up to 3 h before peak flow).Following new estimations of the Jacobian matrix, H, at the and S = 142 mm.The colour scheme is the same as Fig. 7, except for features specific to the pseudo-forecast mode.The black vertical line represents the end of the assimilation period and the start of the forecast period (3 h before the flood peak).The first iterate of the outer loop is in black.All simulations have 5 iterates of the outer loop, however, the algorithm converges after the second iterate in this case, so only the first iterate and the final analysis are shown.
analysis location, the final NS after all iterations of the outer loop is 0.62.The final α was 0.61, suggesting that the algorithm underestimates the rainfall in pseudo-forecast mode when compared to reanalysis mode.The analysis hydrograph is still improved over the background hydrograph, as it reduces the amount of rainfall; however, the reduction is overestimated when only the start of the episode is assimilated.

Results and discussion
This section presents the results of data assimilation applied in 2 modes: reanalysis and pseudo-forecast.The results of the reanalysis mode are discussed in Sect.3.1, followed by the pseudo-forecast mode results in Sect.3.2.In reanalysis mode, results are compared to the background simulation, then to simulations forced with the MFB-corrected rainfall.An analysis of situations in which the algorithm failed to provide an improvement in the discharge forecast concludes Sect.3.1.The pseudo-forecast results start with an analysis of possible sources of error in this mode, followed by the results of assimilation using different parameterisations of the potential storage depth S.

Impact of the rainfall correction
Figure 9 presents NS values in the reanalysis mode compared to the background state for 19 episodes with 7 additional peaks due to separation of multi-peak episodes.The NS values for simulations using uncorrected radar rainfall (the background simulation) are poor and in most cases are not of sufficient quality to reproduce the flood event.Compared to the background state, the NS values of the analysis simulations are improved by an average of 0.75 and are between 0.5 and 1 for a majority of episodes with an average of 0.70.PH values were improved by −0.39 on average (improvements are negative for PH which has an optimal value of 0) and have an average value of 0.14 following assimilation.85 % of episodes show improvement compared to the background state with 15 % showing neutral or negative change following data assimilation.The only degraded episode is that of December 2003; this deterioration is related to the 300 m 3 s −1 upper assimilation limit described in Sect.2.4 and is discussed in greater detail in Sect.3.1.3.Following data assimilation, radar rainfall is of suitable quality for hydrological simulation in most cases.The next section will focus on the comparison of data assimilation to another multiplicative corrector of radar rainfall, the MFB.

Comparison of data assimilation to the MFB correction
A linear regression was performed between MFB and α values for past rainfall episodes as shown in Fig. 10.The two quantities are expected to be related as they both represent corrections of the same rainfall.If errors due to other sources are minimised (parameterisation of the model, measurement of the rain gauges and discharge), the two corrective factors should tend towards the same value.These two quantities are well correlated with a R 2 equal to 0.77.The slope, however, is 1.12, which suggests a systematic underestimation of the rainfall by the MFB correction if α is considered to be the optimal state.The difference between the simulated discharges resulting from the rainfall corrected by the DA procedure and the  In most cases, α provides improved results over the MFB correction.However, some of the improvement in the simulations with α when considering double peaks may be due to an increased time resolution.The MFB was calculated using rainfall over the entire event, whereas the events were separated into single peaks when using α.The MFB is also calculated over a much larger spatial extent than that of the physical basin, leading perhaps to representativeness errors.

Limitations of the assimilation technique
The quality of the December 2003 simulation (Fig. 12a) was degraded following data assimilation when compared to the background state.This is the result of a non-monotonic error in the discharge during the episode, as seen in Fig. 12b.Positive errors in the rising and descending limbs of the hydrograph result in an analysis state with a reduced rainfall.However, the sign of the error in the region near the peak is negative and this part of the hydrograph is not well-represented.To counteract this problem, the upper limit of assimilated observations can be increased to include more observations at the hydrograph peak.The inclusion of these points increases the number of negative errors taken into account by the algorithm and results in an analysis which decreases rainfall less than when discharge observations are limited to less than 300 m 3 s −1 .

Analysis of different sources of uncertainty
In pseudo-forecast mode, the efficiency of the DA algorithm is affected by both a lack of information about the event (representativeness errors) and a poor parameterisation compared to the a posteriori S values (S cal ).Representativeness errors refer to the fact that the start of the event may not be indicative of what comes later.For example, the algorithm would miss the peak region if it were to match observations at the start of the event as closely as possible.Testing a range of β obs values helps to estimate the uncertainty coming from the observations (representativeness), while the comparison of the data assimilation results using different S values gives an idea of the uncertainty resulting from the parameterisation.
To compare the effects of the two sources of uncertainty discussed above, NS and PH values were compared for simulations calculated in pseudo-forecast mode with (i) parameterisation using S Hu2 (β obs = 0.25 m 6 s −2 ) and (ii) different values of the R matrix (β obs = 0.25, 25 and 250 m 6 s −2 ) and S cal .Figure 13a presents a box plot of the change in NS for the four cases and Fig. 13b presents the results for PH.The error in the parameterisation affects the median, as seen by the decreased median for the simulations using S Hu2 , while representativeness errors affect the spread of the results, as seen by the changing width of the distribution for different values of β obs .While β obs serves to limit the influence of observations which do not well represent the rest of the event, it does not bring any new information to the assimilation system.To get a better understanding of what is lost when the event is not fully described, the reanalysis mode can provide an idea of how the information contained in the complete event hydrograph affects the assimilation results.Errors in representativeness are estimated by comparing the difference in NS between the background and analysis simulations in pseudo-forecast mode (S cal , β obs = 0.25) and in reanalysis mode (S cal , β obs = 0.25).The average improvement in NS is 0.35 in pseudo-forecast mode, compared to 0.75 in reanalysis mode.The improvement possible using data assimilation is, thus, cut in half when only the start of the event (until 3 h before the peak) is considered.This process would likely be further complicated if applied in a real-time forecast environment because the peak arrival time would be unknown.
Because representativeness errors are a significant source of uncertainty, β obs was selected as 250 m 6 s −2 for tests using different catchment wetness state indicators to initialise S. As seen in Fig. 13, this value of β obs helps to limit the extent of the change in performance criteria into the negative range.

Results for 3 different soil moisture parameterisations
Figure  NS results may be more positive than the PH results because NS takes into account the assimilation and forecast periods.
In addition, it should be noted that the DA algorithm seeks to reduce the distance between the observed and simulated hydrographs as a whole and not simply at the peak region, thus, it is not expected that DA will always improve peak criteria.
For the NS criterion, 67, 71 and 67 % of episodes were improved by DA using the Claret, Bois Saint Mathieu and Hu2 parameterisations, respectively.For the PH criterion, 67, 62 and 64 % of episodes were improved by DA using the Claret, Bois Saint Mathieu and Hu2 parameterisations, respectively.In general, the three catchment wetness state indicators had similar performances with a slight preference for Claret, which has a higher average Nash value and a tighter PH distribution than the other two indicators.Despite expectations that Hu2 would be the lowest performing catchment wetness state indicator, there is little evidence that modelling errors introduced by this indicator are more important than the uncertainty associated with the two piezometers.
Regressions were performed between the α values and the MFB for each catchment wetness state indicator in addition to S cal .As in Sect.3.1.2,α and MFB are expected to tend towards the same value if uncertainties are minimised.The coefficients of determination, slopes and y-intercepts are presented in Table 4. Contrary to what might be expected, S cal does not have the highest coefficient of determination.This can be in part explained by the random, time-varying nature of radar rainfall and its impact on discharges, which is one of the causes of the representativeness errors mentioned earlier.
In addition to the possible influence of the karst, random errors in the radar rainfall may lead the algorithm to predict a rainfall correction during the start of the rainfall event which does not hold true for the rest of the hydrograph.On the other hand, S cal has a slope of 0.95 and a y-intercept of 0.00 compared to a slope of 0.77 to 0.79 and a y-intercept of 0.07 This should be expected as S cal already contains information about the rainfall gathered through the calibration process.Bois Saint Mathieu, Claret and Hu2 all had similar slopes, which may point to a tendency of the algorithm to underestimate the rainfall correction when initialising the model with measures of the catchment wetness state.All of the regressions had relatively poor coefficients of determination, 0.36 to 0.47, with Claret having the highest value.This is likely due to the random, time-varying nature of the errors in the radar rainfall.The correction calculated using DA will reflect the optimal rainfall multiplier for the start of the rainfall event, which may differ from the MFB correction which is averaged over all event time steps.In the case of S cal , this correction tends toward the MFB correction, but is affected by representativeness errors introduced through random variations in the rainfall during the event.
These results highlight the challenges associated with using a conceptual hydrological model to forecast flood events given the need for model initialisation.The poor quality of the coefficients of determination is an important reminder of the impact of random, time varying error in the rainfall together with uncertainty in the model representation of complex physical processes.Despite the presence of random, time-varying errors, S cal did have improved average NS and PH values, as expected.
Few conclusions can be drawn from the comparison of the different catchment wetness state indicators.The modelled indicator, Hu2 had a performance similar to that of the piezometers.Thus, Hu2 contains information about the catchment wetness state comparable to that of the piezometers and both of the piezometers selected provided adequate information on the catchment wetness state.

Summary and conclusions
A non-sequential Extended Kalman Filter (EKF) was implemented on top of a distributed, event-based, parsimonious rainfall-runoff model.Discharges observed at the catchment outlet were assimilated in order to correct radar rainfall inputs using a multiplier (α) held constant during a given event.The data assimilation (DA) algorithm was effective in both reanalysis and pseudo-forecast modes, despite increased uncertainty due to representativeness and parameterisation errors in the later.Improvements in the model structure might be capable of increasing the efficiency of this DA system, but modelling karstic catchments remains a significant challenge for hydrologists and lies outside the scope of this study, which focuses primarily on the utility of DA for correcting rainfall measured by weather radar.
In reanalysis mode, the DA algorithm is capable of finding an optimal control vector that produces simulations improved over those produced by the mean field bias (MFB) for most episodes given an appropriate parameterisation.These corrections are well correlated with MFB values.
In pseudo-forecast mode, over 60 % of episodes had improved Nash-Sutcliffe efficiency criteria (NS) following data assimilation.Average improvement in the NS was notable, while that of the PH was near 0. These results were subject to representativeness and parameterisation errors which diminished the efficiency of data assimilation compared to the reanalysis mode.Representativeness errors were estimated by comparing the performance of the algorithm in reanalysis and pseudo-forecast modes.Nash-Sutcliffe criteria were improved by an average of 0.35 in pseudo-forecast mode, compared to 0.75 in reanalysis mode, demonstrating that corrections predicted during the start of the event may not be optimal for reproducing event hydrographs.Errors in representativeness may be due to the time-varying nature of the uncertainty in the radar rainfall or model difficulties in representing physical processes in the catchment.To estimate the error resulting from the parameterisation, data assimilation was performed in pseudo-forecast mode with S cal (model initialised using ground and high quality radar rainfall) and then compared to results using S initialised with catchment wetness state indicators.On average, improvements in the NS and PH values of simulations using S cal are nearly double of those initialised using wetness state indicators.It was also seen that the α values from tests using S cal were closer to the MFB than tests initialised using the indicators.Information contained in the model initialisation may help the algorithm to find a correction which reproduces the effect of the MFB.However, regressions between α and MFB values had poor coefficients of determination for all S initialisations due to the representativeness errors which affect the assimilation results in the pseudo-forecast mode.
In both reanalysis and pseudo-forecast modes, errors in simulated discharges occurred due to simplifications of the physical system in the model representation and poor knowledge of the karstic aquifer.Errors also resulted from variations in the rainfall error during the episode, since data assimilation is performed using a constant rainfall correction.This is especially pertinent for the pseudo-forecast mode because only the start of the event is known.A sliding assimilation window or an autoregressive update function may be necessary to improve the analysis quality when the rainfall error varies during the episode.The use of a sliding window to calculate α, with comparisons made to MFB values calculated with the same temporal resolution could help to estimate the efficiency of such a technique.Using a distributed rainfall correction is another possible approach, given the sensitivity of radar measurements to distance (Kahl and Nachtnebel, 2008).However, the updating procedure used in Kahl and Nachtnebel (2008) is limited in that it uses an objective function which does not account for errors in the observations.The tests carried out in the pseudo-forecast mode in this study have shown that the observation error must be accounted for when the event is not completely known.
From a prevision standpoint, testing modelled future rainfall with this algorithm is essential for judging its utility for operational flood forecasting.At the present time, modelled rainfall is not available at a suitable temporal resolution for this region.Further research would also be necessary to adapt this technique for other types of models and floods.This case relates to a conceptual, event-based model used for flash flood events, but physically-based models may prove to be more robust in forecast environments when sufficient data on the watershed is available.As seen in this study, event-based models have the disadvantage of being strongly dependent on the initialisation selected.Floods based on phenomena which take place at a longer timescale may also lead to different results.
In spite of certain limitations of this assimilation system, it may be useful for the correction of radar rainfall following a careful calibration of model parameters.For basins that have available radar rainfall, but scarce or inaccurate ground rainfall measurements, discharge measurements could serve as a replacement for the MFB correction using an appropriate hydrological model and assimilation procedure.

Fig. 1 .
Fig. 1.Visualisation of the Lez catchment and its monitoring network: map of the Lez catchment, the rain gauges used for the measurement of ground rainfall, and the Nîmes weather radar.
to the watershed outlet by what is referred to here as the transfer function.The Lag and Route transfer function (Tramblay E. Harader et al.: Data assimilation for the correction of radar rainfall et al., 2011

Fig. 3 .
Fig. 3. Discharge as a function of α at 3 h before the flood peak.

Fig. 4 .
Fig. 4. Schematic representation of the hydrological model: inputs (blue), model parameters (purple), state variables (dark red) and the background model outputs (Q sim,b in pink).Inputs, parameters, state variables or model outputs can be corrected by DA using observations (Q obs in light blue) and model outputs in order to produce the corrected discharge (Q sim,a in green).

Fig. 5 .
Fig. 5. Schematic representation of the EKF: the background model trajectory (Q sim,b in pink) is corrected using observations (Q obs , blue crosses) to produce the analysis model trajectory (Q a in green) during steps 1 through 4. In step 5, the observation operation is re-linearised in the vicinity of the analysis and steps 2 and 3 are repeated to form an "outer loop".

Fig. 6 .
Fig. 6.The outer loop process.The x-axis represents the value of the control vector and the y-axis is the misfit cost (cost function).The red curve represents the non-quadratic true value of the cost function, while the dotted curves represent successive iterations of the outer loop, each with a new estimate of the Jacobian of H in the vicinity of the previous analysis.

Fig. 7 .
Fig. 7. Reanalysis mode, November 2008: β obs = 0.25 m 6 s −2 and S = 142 mm.The horizontal dashed line is the lower assimilation threshold (2 m 3 s −1 ).Observations are in blue, the background simulation in pink and the analysis simulation in green.Assimilated observations are marked with blue crosses.The hyetogram is on the inverted y-axis: initial rainfall is in dark blue and the corrected rainfall is in light blue with each bar the width of a 1 h time step.This colour scheme is conserved throughout the paper.

Fig. 8 .
Fig. 8. Pseudo-forecast mode, November 2008: β obs = 0.25 m 6 s −2and S = 142 mm.The colour scheme is the same as Fig.7, except for features specific to the pseudo-forecast mode.The black vertical line represents the end of the assimilation period and the start of the forecast period (3 h before the flood peak).The first iterate of the outer loop is in black.All simulations have 5 iterates of the outer loop, however, the algorithm converges after the second iterate in this case, so only the first iterate and the final analysis are shown.

Fig. 9 .
Fig. 9. Comparison of background NS values with NS values following data assimilation (analysis).The x-axis contains the episode label in the format mYYpp, where m is the first letter of the month (j is January and m is May), YY is the year and pp is the peak number for the 2nd and greater peaks.

Fig. 10 .
Fig. 10.Regression of α versus MFB.y = x is drawn in red and the regression in blue.

Fig. 11 .
Fig. 11.Improvements in simulation quality indicators for the 19 rainfall episodes.The dark blue bars represent NS α − NS MFB .The light blue bars are the difference in the normalised peak flow criteria, PH MFB − PH α .

Fig. 12 .Fig. 13 .
Fig. 12. Reanalysis mode, December 2003: (a) discharges: observations are in blue, the background simulation in pink and the analysis simulation in green; (b) the error in the simulated discharge, Q background − Q observations (red).

Table 1 .
Rainfall events occurring over the Lez catchment from 1997-2008.The date, mean field bias (MFB) and peak discharge (Q peak ) are shown.

Table 2 .
The S -catchment wetness state indicator relationship.M is the slope of the linear regression between S cal and the wetness state indicator, b is the y-intercept, and R 2 is the coefficient of determination for this regression.% change refers to the average difference between S reg and S reg using the validation period regression.σ is the standard deviation of this difference.

Table 3 .
S reg estimated from physical indicators of the catchment wetness state using the linear regressions presented in Table2.Dashes indicate missing values.
14 presents box plots of the improvements in the NS and PH values for the three different S parameterisations.β obs is selected as 250 m 6 s −2 .Bois Saint Mathieu and Claret both have an increased median NS improvement compared to Hu2.The spreads of Bois Saint Mathieu and Claret improvements are similar.For the PH criterion, the medians of each of the three catchment wetness state indicators are similar, though Claret has the narrowest spread, but also several negative outliers.The NS was improved by an average of 0.23, 0.31 and 0.16 for Bois Saint Mathieu, Claret and Hu2, respectively, compared to 0.40 for S cal .The PH was improved by an average of 0.07, 0.04 and 0.07 for Bois Saint Mathieu, Claret and Hu2, respectively, compared to 0.14 for S cal .The Hydrol.Earth Syst.Sci., 16, 4247-4264, 2012 www.hydrol-earth-syst-sci.net/16/4247/2012/

Table 4 .
α-MFB regression for catchment wetness state indicators.Notation follows that of Table 2.
cal are thus much closer to the MFB values than those of the indicators if we were to consider the regression alone.