STATE UPDATING OF A DISTRIBUTED HYDROLOGICAL MODEL WITH ENSEMBLE KALMAN FILTERING: EFFECTS OF UPDATING FREQUENCY AND OBSERVATION NETWORK DENSITY ON FORECAST ACCURACY

This paper presents a study on the optimal setup for discharge assimilation within a spatially distributed hydrological model. The Ensemble Kalman filter (EnKF) is employed to update the grid-based distributed states of such an hourly spatially distributed version of the HBV-96 model. By using a physically based model for the routing, the time delay and attenuation are modelled more realistically. The discharge and states at a given time step are assumed to be dependent on the previous time step only (Markov property). Synthetic and real world experiments are carried out for the Upper Ourthe (1600 km 2 ), a relatively quickly responding catchment in the Belgian Ardennes. We assess the impact on the forecasted discharge of (1) various sets of the spatially distributed discharge gauges and (2) the filtering frequency. The results show that the hydrological forecast at the catchment outlet is improved by assimilating interior gauges. This augmentation of the observation vector improves the forecast more than increasing the updating frequency. In terms of the model states, the EnKF procedure is found to mainly change the pdfs of the two routing model storages, even when the uncertainty in the discharge simulations is smaller than the defined observation uncertainty.


Introduction
Accurate and reliable hydrological forecasts have been a challenge in applied hydrology for decades.Better forecasts can be obtained through data assimilation (DA) by merging observations with model simulations (Reichle, 2008).This approach basically updates the model states with externally measured variables (Pauwels and De Lannoy, 2006) to obtain correct initial conditions for the next time step.Currently, most operational hydrological forecasting systems employ lumped hydrological models (with deterministic or manual state updating), but there is a clear tendency to move towards distributed models combined with hydrological ensemble forecasts, (e.g.Koren et al., 2004;Cole and Moore, 2009;Weerts et al., 2012).The main advantage of spatially distributed models is the possibility to force them with spatially measured data, which nowadays become more readily available due to rapid developments in telemetry.Distributed model states also resemble the real world observations (e.g., groundwater levels, soil moisture, discharge) at the interior of the catchment more closely than lumped states over the whole catchment.Another advantage of applying distributed models is the ability to simulate and predict hydrological variables at interior locations within the catchment.Techniques on how to perform ensemble data assimilation using these models in real-time settings should be developed and tested (Weerts et al., 2012;Liu et al., 2012).
Data assimilation methods used in hydrology can be divided into two classes: (1) sequential and (2) variational (e.g.Liu and Gupta, 2007).Sequential methods are mostly employed for state updating in hydrological models by assimilating observations for each time step when they become O. Rakovec et al.: State updating of a distributed hydrological model with EnKF available.Its impact depends on the uncertainties in both the observations and model states.Variational methods rather minimise a cost function over a simulation time window.At the beginning, a first-guess model is constructed, which is afterwards updated by creating an adjoint model which propagates backwards in time and incorporates the mismatch between the model and observations (Liu and Gupta, 2007).
A popular method often used in both meteorological and hydrological forecasting is the Ensemble Kalman Filter (EnKF) (Evensen, 1994(Evensen, , 2003(Evensen, , 2009)).This sequential data assimilation method is an extension of the classical Kalman filter (KF), which was originally developed for linear systems (Kalman, 1960).The EnKF propagates an ensemble of model realisations through time and estimates the error covariance matrix from the ensemble statistics.The advantage of the EnKF to other data assimilation methods is computational efficiency, and easy and straightforward implementation within a data assimilation procedure for both lumped and distributed models (Pauwels and De Lannoy, 2009).
Discharge measurements are the most widely used in situ hydrological observations for model updating, since they reflect the local catchment wetness conditions and are often available at high temporal resolution (Pauwels and De Lannoy, 2006;Teuling et al., 2010), which is necessary for operational hydrological forecasting.Comprehensive analyses of discharge data assimilation into spatially lumped hydrological models was carried out by e.g., Pauwels and De Lannoy (2006); Weerts and El Serafy (2006); Pauwels and De Lannoy (2009).Among others, Aubert et al. (2003) and Lee et al. (2011) assimilated, in addition to discharges, in situ soil moisture measurements.Lee et al. (2012) analysed the sensitivity of variational data assimilation methods for multiple spatiotemporal adjustment scales, namely assessing (1) different spatial distributions of model states and (2) temporal resolution of biases in precipitation and potential evaporation.Furthermore, the water level is an example of an in situ state variable, which can also be assimilated into a hydrodynamical model (e.g.Madsen and Skotner, 2005;Neal et al., 2007;Weerts et al., 2010).Additionally, temperature observations play an important role in real-time operational forecasting systems, especially in regions with significant snow melt (e.g.Verbunt et al., 2007;Sene, 2008).
Although the assimilation of remotely sensed data into operational hydrological models can improve model performance, this task is complicated, because remotely sensed data usually have a higher uncertainty than other in situ measurements (Pauwels and De Lannoy, 2009).Nevertheless, several studies focus on assimilation of remotely sensed data into hydrological models, e.g., soil moisture (Pauwels et al., 2001(Pauwels et al., , 2002;;Moradkhani, 2008) and snow (Slater and Clark, 2006;Thirel et al., 2011).Remotely sensed water levels were assimilated into a hydraulic model by Giustarini et al. (2011).
So far, few studies have reported on data assimilation within spatially distributed hydrological models, both from a scientific and from an operational perspective.Clark et al. (2008) assimilated discharges using the adjusted EnKF for a real world case only.However, their approach was not tested and evaluated in a synthetic experiment, which would help to understand the behaviour of the simulated and updated probability density function (pdf) of the model states.Another EnKF study, by Blöschl et al. (2008), applied discharge assimilation into a real-time grid-based operational flash flood forecasting model.As Blöschl et al. (2008) did not elaborate on the importance of individual discharge gauges in the interior of the catchment, this remains an interesting question to address.However, more recently the positive effect of interior discharge gauges on hydrological forecast was described by Lee et al. (2011), who assimilated both discharge and/or in situ soil moisture data using the variational method.
From an operational point of view, it is also desirable to optimise the telemetry system, which delivers the observations to the forecaster.One of the interesting aspects is the frequency at which the observations become available.In most of the currently operational large-scale flood forecasting systems, it is typical to separate the forecast process into a state (i.e., carry-over or update) run, which is normally run once a day, and forecast runs which are run more frequently.The frequency of the state run is very much dictated by the frequency at which observed data becomes available (i.e., frequency of polling of the telemetry system).Hence, it is interesting to test the optimal updating frequency at which the hydrological data are assimilated into the forecasting model to obtain the most accurate forecast.
The objective of this study is to analyse the sensitivity of the DA procedure to the number and the locations of discharge gauges, which are assimilated into a grid-based distributed operational hydrological forecasting model using the EnKF.The optimal updating frequency will be addressed as well.
The focus here lies mainly on the state run part of the forecast chain and, therefore, we leave out the meteorological forecasts and error modelling of the meteorological forecast, which is a research topic on its own (e.g.Germann et al., 2009).Using this approach, we disregard the large uncertainties in the meteorological forecasts and we can focus entirely on the uncertainty in the initial model states.

Catchment description and data availability
The hydrological simulations are carried out for the Upper Ourthe catchment, upstream of Tabreux (Fig. 1b), which drains an area of about 1600 km 2 (Berne et al., 2005).This catchment forms a tributary of the Meuse River basin, originating in the hilly parts of the Belgian Ardennes.The climatic conditions can be classified as rain-fed with irregular snow in winter and the runoff regime is highly variable with a Location of discharge gauges is indicated in Fig. 1b.
low summer discharges and high winter discharges (Leander et al., 2005).Relatively shallow soils in combination with significant elevation differences result in a quickly responding catchment.As such, the whole region represents a significant flood risk to the Netherlands (de Wit, 2008).Table 1 presents the catchment response time of the Upper Ourthe, which is defined as the time between the event-based catchment averaged rainfall centroid and the discharge peak for individual discharge gauges.The catchment response time is about 30 h at Tabreux and about 11 h at the two most upstream gauges (Mabompre and Ortho).This indicates that it takes about 20 h of travel time within the main channel between the two upstream gauges and the catchment outlet.
Hourly precipitation data are available from 42 automatic rain gauges situated within the Belgian Ardennes region, from which 10 are located inside the Upper Ourthe catchment (Hazenberg et al., 2011).Discharge is measured at six different points at an hourly resolution.Next to that, temperature is obtained from the Saint Hubert meteorological station (Fig. 1b).The long-term mean monthly values of potential evapotranspiration are assumed identical to those of the operational lumped HBV-96 model, derived from the St. Mihiel station in North-Eastern France.

Generation of a spatially distributed precipitation ensemble
A precipitation ensemble, a finite and discrete number of spatial realisations over time, represents the uncertainty associated with temporal as well as spatial variation of precipitation.This probabilistic input enables hydrologists to evaluate more critically hydrological simulations and/or predictions.
In the current paper, we employed a time-dependent multivariate spatial conditional simulation method (Rakovec et al., 2012), which is further made conditional on preceding simulations.This method identifies, at an hourly time step, temporally coherent errors in spatial precipitation fields that are plausible from a hydro-meteorological perspective.Neglecting this temporal aspect could lead to underestimation of the overall uncertainty in the precipitation ensemble.The theory of conditional (sequential Gaussian) simulations is fully explained by Goovaerts (1997).For a detailed description of this time-dependent simulation method using the gstat R package (Pebesma, 2004;Rossiter, 2007;R Development Core Team, 2011) we refer to Rakovec et al. (2012), who carried out a synthetic experiment and analysed three real rainfall events during winter 2002/2003 within the Belgian Ardennes.Altogether, 27 rain gauges were used to simulate 64 ensemble members over a 100 km × 100 km domain with a 10 km × 10 km raster resolution (see Fig. 1b).In the current study, we made the simulation of each realisation conditional on 3 h of the corresponding previously simulated realisation, which is a recommendation following from the results obtained by Rakovec et al. (2012).

Hydrological model
Currently, a spatially lumped HBV-96 model (Lindström et al., 1997) is used operationally by the Dutch authorities for flood forecasting of the Meuse River basin at and downstream of Sint Pieter at an hourly time step (Driessen et al., 2010).The Meuse River basin upstream of Sint Pieter (∼ 21 000 km 2 ), the entrance point into the Netherlands, is conceptualised into 15 lumped HBV-96 sub-catchments of which the Upper Ourthe is one (Fig. 1a).The models for the Meuse River sub-catchments were calibrated for the period 1970-1984 and validated for the period 1985-1996 at a daily time step by Booij (2002).The precipitation inputs for the original HBV-96 models were derived from 39 rain gauges, of which only one station was located within the Upper Ourthe catchment.For operational purposes, the hourly HBV-96 models were derived and re-calibrated based on the work of Booij (2002) by van Deursen (2004).A comparison of daily and hourly model versions was carried out by Weerts (2007).Currently, the operational forecasts derived with the lumped hourly HBV-96 models are used as lateral inflows into a 1D-hydrodynamic model of the Meuse River.This lumped model does not employ sequential state updating, but is updated with discharge observations by means of an automated auto-regressive moving average error correction method (Broersen and Weerts, 2005).
For this study, we have developed a grid-based spatially distributed HBV-96 based model within PCRaster, a software environment for constructing iterative spatiotemporal environmental models (Karssenberg et al., 2009;PCraster, 2012).Such a grid-based approach is a popular concept in applied hydrology (e.g.Koren et al., 2004;Blöschl et al., 2008;Cole and Moore, 2009;Thielen et al., 2009).For each 1 km × 1 km grid cell of the Upper Ourthe the HBV-96 model was implemented and is used as a benchmark case.Both model parameters and structure are taken identical to the lumped HBV-96 version used operationally except for the discharge routing, for which a kinematic wave model (Chow et al., 1988;PCraster, 2012) is used.This is only the first step forward by taking full advantage of distributed models, which would allow us to define spatially variable model parameters.However, that is beyond the scope of the current study.
The structure of the model used in this study is shown in Fig. 2. For each grid cell the model considers four model states: (1) snow (SN), (2) soil moisture (SM), (3) upper zone storage (UZ) and (4) lower zone storage (LZ).The dynamics of the model states are governed by the following model fluxes: rainfall, snowfall, snowmelt, actual evaporation, seepage, capillary rise, direct runoff, percolation, quick flow and base flow.The latter two fluxes force the kinematic wave model.This routing scheme calculates the overland flow using two additional model states, the water level (H) and discharge accumulation over the drainage network (Q).In this study, we use a very similar routing setup as the one applied within the distributed hydrological CQ-flow model (Schellekens, 2006).The main drainage network is obtained using the 8-direction steepest descent algorithm based on a digital elevation model with a grid resolution of 1 km × 1 km.Afterwards, the catchment is split between the channel and non-channel grid cells.The channel network is defined for the cells with Strahler stream order (Strahler, 1964) greater than 3.As such, only the major tributaries of the Upper Ourthe are identified, which corresponds well to the channel network derived from a topographic map (for comparison, see Figs. 1b and 4).A channel width is assigned based on a field survey and Google maps (2011) according to the Strahler stream order number, as shown in Table 2. Google orthophoto maps (satellite images) were used for rough visual estimates to obtain information on the channel widths of corresponding Strahler stream order numbers.For the non-channel cells, the channel width remains equal to the grid cell width, because the water is routed from these cells by means of sheet flow on top of the whole grid cell.By making the width of the nonchannel cells very large, we are able to decrease the hydraulic gradient of the water in the "channel" of these non-channel cells and, therefore, increase the response time of these rather slow responding cells.The channel width is then used to derive the water level, which defines, together with the local topography gradient, the gravity force, which is a driver for the river flow.Channel roughness coefficients were estimated by a sensitivity analysis.
To demonstrate that the implemented spatially distributed version is able to produce reasonable hydrological discharge simulations using the aforementioned precipitation ensemble generator, Fig. 3 shows the observed discharge and both the spatially lumped and distributed HBV-96 discharge simulations at Tabreux.A clear consistency can be observed between both HBV-96 simulations for the 5 month period (15 August 2002-15 January 2003), which will be further used in this study.The Nash-Sutcliffe (NS) model efficiency coefficient (Nash and Sutcliffe, 1970) between both HBV-96 simulations is about 0.99, which gives strong evidence of model similarity.Moreover, the NSs between the observed and the ensemble of discharges using the grid-based version of HBV-96 are between 0.92 and 0.96 and the rootmean-square error (RMSE) ranges between 9 and 12 m 3 s −1 (Fig. 3).The NS between the observed and the simulated discharge using the lumped HBV-96 is about 0.96 as well.For completeness, additional statistics are summarised in Fig. 3.Even though the grid-based HBV-96 model was not recalibrated, the model performs very well at the scale of the Upper Ourthe catchment.

Ensemble Kalman filter
The Ensemble Kalman filter (Evensen, 2003(Evensen, , 2009;;Weerts and El Serafy, 2006) is a recursive Bayesian estimation method, which estimates the true probability density function of the model states conditioned on observations.Let us denote a dynamic state space system as: where x k is a state vector at time k, f is an operator expressing the model state transition from time step k − 1 to k in response to the model input u k−1 and time-invariant model parameters θ.So f is, in fact, the hydrological model.ω k stands for system noise, normally distributed with zero mean and covariance S.This additive system noise incorporates the overall uncertainties in model structure, parameters and model inputs.One can expect some spatial patterns of model errors to be found in the covariance matrix S.However, quantification of S for highly nonlinear hydrological systems is a complicated task and, therefore, we keep it time-invariant.
The observation process is governed by Eq. ( 2): where y k is an observation vector derived from the model state x k and the model parameters through the H operator (in our case the kinematic wave model).ν k stands for additive observational noise, normally distributed with zero mean and covariance R k .For independent measurement errors between the observations in vector y k , we can assume R k to be a diagonal matrix.As such, this simplification does not consider any dependency in model simulations for observation points located close to each other.The idea of recursive Bayesian estimation is to construct a conditional density p for an ensemble of the state x k given all available information up to and including the step k: p(x k |Y k ), where Y k = (y 1 , y 2 , . . ., y k ).This can be obtained using the Bayesian rules in two steps: forecast and update.
After the update of model states at time k−1, the hydrological model is used to forecast model states at time k (Eq.1).The grid-based model states form a matrix, which consists of N state vectors x k corresponding to N ensemble members: where SN, SM, UZ, LZ, H and Q are the HBV-96 model states (Sect.2.3), m gives the number of grid cells and T is the transpose operator.The ensemble mean is used to derive the model error for each ensemble member: The ensemble estimated model covariance matrix P k is defined as When observations become available, the model states are updated as follows: where X + k is the new updated (posterior) model state matrix and X − k is the forecasted (prior) model state matrix.K k is the Kalman gain, a weighting factor of the errors in model and observations: where P k H T k is approximated by the forecasted covariance between the model states and the forecasted discharge, and k is approximated by the variance of forecasted discharge (Houtekamer and Mitchell, 2001): where In previous published papers on EnKF (e.g.Weerts and El Serafy, 2006;Pauwels and De Lannoy, 2009), a time delay issue was noted due to the use of the unit-hydrograph, where the discharge at time k depends on several previous calculated discharges (for instance at k −1, k −2, k −3, k −4, etc., for hourly models often up to 24 h or more in the case of the Upper Ourthe).By using a physically based model for the routing (Sect.2.3), the time delay and attenuation are modelled more realistically and the discharge and states x k can be assumed to depend only on the states x k−1 (Markov property).The time delay is, thus, explicitly taken into account in the model.

Synthetic experiment
The approach of the synthetic experiment is similar to the approach used by Weerts and El Serafy (2006), with the main difference being that we employ a realistic stochastic representation of the spatially distributed precipitation (see Rakovec et al., 2012).We limit the analysis of the experiment to input uncertainty only.The main reason for this is that we want to fully understand and investigate the filter process using a distributed hydrological model and realistic precipitation fields.We believe this is already challenging enough without initial state and model parameter/structural uncertainty, which we leave for future work.
The synthetic experiment has been carried out to examine the ability of the EnKF to update the grid-based hydrological model states via assimilation of the spatially measured discharge.The EnKF procedure was applied to a 5month period from 15 August 2002 to 15 January 2003, which includes a dry and a wet period.For reasons of clarity of the experimental setup, we did not update model parameters, which were kept constant.In an open loop simulation, i.e., without data assimilation, the model was initially forced with uncertain precipitation inputs with a simulation memory of 3 h (64 ensemble members) derived using time-dependent multivariate spatial conditional simulation (see Sect. 2.2 and Rakovec et al., 2012) and observed deterministic potential evapotranspiration and temperature.This produced an ensemble of simulated discharges from which one complete realisation was randomly selected as the true discharge (Q true,k ).To introduce discharge observation uncertainty, a normally distributed error ν k (Eq.2) with heteroscedastic variance was added to Q true,k to obtain a synthetic perturbed observation Q pobs,k : The discharge measurement error was defined similar to Thiemann et al. (2001) and Weerts and El Serafy (2006).
The ensemble size of 64 members corresponds to other studies as well, e.g., Pauwels and De Lannoy (2009).The model errors (S, Eq. 1) were obtained by perturbing the model states indirectly by uncertain precipitation input.Additional direct perturbation of the model storages, which was applied by e.g., Clark et al. (2008) was not considered, because this was considered to be beyond the scope of this study.That would make our example, which focuses on the input uncertainty only, even more complicated.
In the synthetic experiment, we assimilated in total 5 schemes of perturbed discharge observations, expressed by the vectors y k , as shown in Fig. 4. The first case A is identical to a lumped model for the Upper Ourthe, where only the most downstream observation is available.The second case B considers only the two most upstream discharge gauges.The third case C includes three additional observations upstream to case A. The fourth case D contains all six discharge gauges and the fifth case E includes an additional 12 imaginary gauges to the fourth case.
Moreover, the effect of the updating frequency, i.e., how often the observations become available and how often they are assimilated into the model, was analysed for updating frequencies of 24, 12 and 6 h.
The performance of the data assimilation procedure regarding discharge forecasting was then evaluated using the root-mean-square error (RMSE lt ): where lt stands for lead time (lt = 1 h, 2 h, . . ., 48 h) and Q for is a forecasted discharge vector of length J .M is the number of hydrological forecasts, which were issued over the 5month period.To allow comparison of the different updating frequencies between each other, the hydrological forecast was issued every 6 h, i.e., 4 times a day.

Real world experiment
In the real world experiment, we applied the same model forcings as in the synthetic experiment (Sect.2.5.1).The difference with the synthetic experiment was the assimilation of the real discharge observations (Q obs,k ), which were perturbed by a normally distributed observation error with a variance of (0.1 Q obs,k ) 2 (Weerts and El Serafy, 2006;Clark et al., 2008).Analogous to the synthetic experiment, Q obs was used in Eq. ( 14) instead of Q true to calculate the RMSE.The size of the observation vector y k was limited to cases A, B and D (Fig. 4).

Model performance regarding discharge forecast
The long-term RMSE (Eq.14) between the synthetic observed discharge and the forecasted discharge of the 64 ensemble members at Tabreux (Fig. 1b) for the three updating frequencies of 24, 12 and 6 h is shown in Fig. 5.The forecasted discharge without data assimilation gives a constant RMSE of about 5 m 3 s −1 , which corresponds to 16 % error with respect to the mean simulation over the 5-month period.
For the forecasted discharge with data assimilation, there is a reduction in RMSE for all discharge observation vectors, however, the magnitudes differ as follows: (1) for the updating frequency of 24 h, the benchmark case A performs worst of all 5 cases, although the differences between case A and B are marginal, given the small difference in their RMSE values as well as the high NS value in the base-case model simulation.Moreover, there is a gradual decrease in RMSE for the cases with a large number of assimilated gauges.
(2) For the updating frequency of 12 h, there is even a further reduction in RMSE for all the observation vectors.The largest reduction is achieved for case E, in which the RMSE at the lead time of 1 h is 1.4 m 3 s −1 (5 % of the mean observed discharge).Additionally, case A (one gauge at the outlet) outperforms case B (two interior gauges only) for lead times up to about 20 h, which is in line with the channel travel time from the most upstream gauges to the outlet.
(3) For the updating frequency of 6 h, there is not a pronounced improvement in RMSE.This can be expected, because within the 6 h between updating moments, hardly any rainfall is transformed into discharge, even at the most upstream gauges, as is shown in Table 1.

State updating
A further logical step in the analysis of the results is to have a look at the ability of the DA procedure to correctly update the model states.In other words, we wanted to check if the setup of the EnKF can identify the pdf around the true model states.However, we recall that there is not a single configuration of model states yielding one discharge value due to the fact that our system is spatially distributed.
We investigate the effect of the observation vector at the updating frequency of 24 h.We selected 2 locations (location 1 and 6 in Fig. 1b) within the catchment domain for which the simulated and updated model states are presented in Figs. 6 and 7. ] q q q q q q q q q q q q q q q q 24h 0 10 20 30 40 Lead time [h] q q q q q q q q q q q q q q q q 12h 0 10 20 30 40 q q q q q q q q q q q q q q q q 06h q No update A B C D E model states for the five discharge observation vectors (cases A, B, C, D and E, see Fig. 4).The true model states are indicated by asterisks.Recall that the true state for the synthetic experiment was a randomly selected sample.Note that the snow model states are not shown because there was no snow simulated during the 5-month period.We have chosen 3 November and 2 January, because both dates occur shortly before a discharge peak, although the wetness conditions of the catchment are different.The first and smaller peak is observed when the model storages are rather dry.The second and larger peak occurs after an extensive rainy period, when the model states became fully saturated.Figure 6 indicates that at the catchment outlet there is hardly any difference between the forecasted and updated model states in soil moisture (SM), upper zone storage (UZ) and lower zone storage (LZ) for all discharge observation vectors and for both dry and wet conditions.However, note that both the forecasted and updated pdfs of SM, UZ and LZ tend to have more accurate peaks around the "true" values for a larger number of assimilated discharge gauges.This means that, even though there is no clear difference between the forecasted and updated pdfs at one time instant, its accumulation over time makes it visible in the higher kurtosis.Therefore, it makes sense to update those rather insensitive model states.Furthermore, the EnKF is well able to identify the "truth" in two routing storages, the water level (H) and the water storage in the channel (Q) on 3 November as well as on 2 January.There is a clear shift of the updated histogram centroid towards the true value for all discharge observation vectors, except for case B, which remains unchanged.This is caused by the fact that case B (Fig. 4) consists of two discharge gauges, which are located far away from the catchment outlet.Furthermore, it can be seen that the EnKF is well able to identify the two routing states even if the prescribed discharge observation error bands are larger than the ensemble spread of the forecasted discharge.
At the interior point (Fig. 7), similar to the catchment outlet (Fig. 6), there is no pronounced update in forecasted and updated histograms for soil moisture (SM), upper zone storage (UZ) and lower zone storage (LZ).The two routing states are again easier to identify.However, the ability of the EnKF to identify the true H and Q depends on the location of stations in the discharge observation vector.For the cases B, D, and E, which contain at least one gauge situated close to location 6, the updated histogram of H and Q approaches the true state and also its shape becomes more peaked.On the other hand, for the cases A and C, which do not include gauges close to location 6, no changes in H and Q histograms occur.

Model performance regarding discharge forecast
The long-term RMSE for the real world experiment is shown in Fig. 8. Similar to the synthetic experiment, all three discharge observation vectors assimilated into the model improve the forecasted discharge for lead times up to 48 h, except for case B, which slightly deteriorates the forecast performance in comparison with the forecasts without discharge assimilation for lead times longer than 30 h.The best results, meaning the lowest RMSE, are achieved by assimilating all six gauges (case D) for all updating frequencies, although for longer lead times it approaches the benchmark case A (the outlet gauge only).The largest reduction was achieved by case D, for which the RMSE at the lead time of 1 h was about 6 m 3 s −1 (20 % of the mean observed discharge) for an updating frequency of 6 h.Similar to the synthetic experiment, case A (one gauge at the outlet) outperforms case B (two interior gauges only) for all updating frequencies in the real world experiment.Moreover, we can observe for case B a rather constant RMSE during the first 20 h.This surprisingly steady RMSE may be explained by the assimilation effect of the most upstream gauges (locations 5 and 6), for which it takes about 20 h to reach the outlet (location 1).Although an increase in updating frequency from 24 h to 12 h improves the RMSE, a further increase in the updating frequency from 12 h to 6 h yields a more or less equal RMSE, which corresponds to the synthetic experiment.
The short-term RMSE for an individual major flood peak, which was observed at the beginning of January 2003, is shown in Fig. 9.Because of the rather short period used in this analysis, the shapes of the RMSE are not smoothed out and the forecasted RMSE without updating is not constant over time either.The best forecast improvements are again ] q q q q q q q q q q q q q q q q 24h 0 10 20 30 40 Lead time [h] q q q q q q q q q q q q q q q q 12h 0 10 20 30 40 q q q q q q q q q q q q q q q q 06h q No update A B D  ] q q q q q q q q q q q q q q q q 24h 0 10 20 30 40 Lead time [h] q q q q q q q q q q q q q q q q 12h 0 10 20 30 40 q q q q q q q q q q q q q q q q 06h q No update A B D achieved by assimilating all six discharge gauges (case D) for all updating frequencies for lead times up to about 15-20 h.For longer lead times, case B (only two upstream gauges) gives very similar RMSE to case D, because the added value of the more downstream gauges (1-4 in Fig. 1b) is filtered out after about 20 h, as shown in Table 1.It is worth mentioning that case B outperforms case A for lead times from 5 h to 20 h, which is not observed using the long-term statistics (Fig. 8).This is caused by the spatial properties of this major flood peak, during which the importance of the upstream gauges is clearly shown, however, completely averaged out in the long-term statistics.

State updating
Like in the synthetic experiment, hardly any change between the forecasted and the updated histograms is observed for soil moisture, upper zone storage and lower zone storage (Figs. 10 and 11), but visible changes can be seen in the routing storages, water level and discharge.For the discharge observation vectors, which contain at least one gauge in the vicinity of the location of the state observation, there is a shift of the centroid of the histograms for discharge and the corresponding water level towards the uncertain discharge observation constrained by the error bars.

Validation
In order to validate the results of the real world experiments, the Nash-Sutcliffe model efficiency (NS) (Nash and Sutcliffe, 1970) is calculated at those stream gauge locations where the discharge data are not assimilated.The median NS of the 64 ensemble updated discharge simulations (Cases A and B) as well as the open loop simulations (without assimilation) are shown in Table 3 for the updating frequency of 24 h.The validation results indicate that the assimilated simulations give performance of the same order of magni- tude as the open loop simulations.Small improvements can still be observed when the assimilated observations are upstream of the validation stations (Case B).However, when the assimilated observation is downstream of the validation gauges (Case A), the interior simulations can slightly deteriorate, but this is observed only at one gauge (number 5) out of five gauges.Nevertheless, this seems to be rather acceptable when we consider the large distance between the assimilation gauge (number 1) and validation gauge (number 5) and the fact that the two most upstream gauges are located at two parallel river branches (see Fig. 1a).

Discussion
The  delay between model states and discharge, as would be needed in spatially lumped models using the retrospective EnKF (Pauwels and De Lannoy, 2006).Another advantage of a grid-based approach over a lumped one is that the spatially distributed discharge observations can be easily incorporated into the model states and make the forecast more accurate for longer lead times.
A novel approach of this study was the application of time-dependent multivariate spatial conditional simulations (Goovaerts, 1997;Pebesma, 2004;Rakovec et al., 2012) of hourly rain gauge observations, used to force the hydrological model in hindcasting mode.As demonstrated by Rakovec et al. (2012), this multivariate approach satisfies for each precipitation realisation the requirement to have a coherent temporal evolution (required within the DA procedure), unlike the time-independent univariate simulations.Using this precipitation ensemble generator, we can achieve the goal that the corresponding simulated spatially distributed model states inherit the temporal aspect of the rainfall fields.As an alternative to rain gauge observations, the precipitation ensemble could be derived from radar rainfall estimates from the C-band radar located in the catchment (Hazenberg et al., 2011), which is a possible topic of further studies.
This study also provided a closer look at the pdfs of the forecasted and updated model states during two hydrologically different situations, while the majority of hydrological DA papers on state updating focus only on the forecasted and analysed discharge and do not address the importance of individual model states.In this study, mainly the pdfs of the two routing model storages were affected by the Kalman filter update, while the other model states (SM, UZ, LZ) were found to be less sensitive to the EnKF procedure.This is because the current formulation of the EnKF (see Eq. 4) does not explicitly consider the strong correlation between soil moisture states in the immediate past and streamflow at the time of forecast.Therefore, it may be difficult to build a covariance matrix among the water balance model states (i.e., SM, UZ, LZ) via assimilating discharge observations.Based on our results, we can state that, given a measured discharge downstream, it is difficult to adjust (and justify) the soil moisture upstream (in a spatially distributed coherent manner) using an EnKF.Other filters like the Ensemble Kalman Smoother (EnKS), which calculate the analysis from several previous time steps (Evensen and Van Leeuwen, 2000), may result in better adjustment of the spatially distributed soil moisture states, which may improve forecasts for even longer lead times.However, with a larger number of assimilated discharge gauges, both the forecasted and updated pdfs of SM, UZ and LZ had more accurate peaks around their true values.Therefore, it makes sense to update those rather insensitive model states.The reason for this behaviour might be the limited model structure, which is similar to other PCRaster operational hydrological applications like the LISFLOOD model (de Roo et al., 2000;Salamon and Feyen, 2009), where the individual neighbouring model cells are not connected by means of interflow and regional groundwater flow, but only drained by some sort of sheet flow via the routing states.This means that the SM, UZ and LZ model states are only controlled by the spatial variation of rainfall.
Finally, it is interesting to note that, unlike Clark et al. (2008), we were able to improve hydrological forecasts using the standard EnKF implementation in both synthetic and real world experiments compared to open loop simulations.

Conclusions and recommendations
We analysed the sensitivity of the data assimilation procedure to the updating frequency, the number and the locations of interior discharge gauges, which were assimilated into a grid-based distributed hydrological forecasting model using the EnKF.By using a physically based model for the routing, the time delay and attenuation are modelled more realistically than when using a unit-hydrograph approach for of the routing.As a consequence the discharge and states at time k can be assumed to depend only on the states at time k − 1 (Markov property).The validation station of this study is Tabreux, which is the Upper Ourthe catchment outlet, Belgium.
In the synthetic experiment we showed that the hydrological forecast at the catchment outlet is improved (in terms of the forecasted root-mean-square error RMSE) by assimilating more interior gauges.This is logical, because all other discharge observations contain information from upstream to improve the posterior forecast.In addition, the EnKF proce-dure is mainly changing the pdfs of the two routing model storages, even when the uncertainty in the discharge simulations is smaller than the defined observation uncertainty.This is because the current formulation of the EnKF does not explicitly consider the strong correlation between soil moisture states in the immediate past and streamflow at the time of forecast.Moreover, with an increasing number of discharge observations, the centroid of the updated histograms within the observation error bounds was approaching the true value more closely and with smaller variance than for the less dense discharge observation networks.
In the real world experiment, the best results in terms of the RMSE were achieved using all observations, which includes all six discharge gauges.Given the travel time of the catchment, an updating frequency of 12 h seems to be the most appropriate.Additionally, similar to the synthetic example, only the two routing model storages showed some sensitivity to the EnKF procedure in terms of the forecasted and updated histograms.We can conclude that the hydrological forecast at Tabreux can be improved by assimilating more upstream gauges using the EnKF data assimilation procedure.This augmentation of the observation vector improves the forecast more than increasing the updating frequency.
For operational use, we recommend the implementation of additional upstream gauges into the observation system, which would enable an increase of the updating frequency and more accurate forecasts, if the polling frequency allows for it.Another recommendation for future research is to have a closer look at alternative model structures (including recalibration of spatially distributed model parameters) and their effect on the sensitivity of individual model states within the EnKF.The main limitation of the current model structure is that there is no flux between neighbouring cells except for the two routing model states.Alternatively, hydrological forecasts can be improved by applying other Kalman-type methods, e.g., the Ensemble Kalman Smoother (Evensen and Van Leeuwen, 2000), which calculates the analysis from several previous time steps up to the time of forecast, instead of mapping the instantaneous covariance between states and discharge (Clark et al., 2008), as shown for the standard EnKF in this study.Finally, additional in situ observations can be considered to be assimilated into the spatially distributed model states, e.g., soil moisture and/or groundwater levels.The latter are believed to resemble point-wise the actual regional water storage more closely than soil moisture observations.the routing scheme in PCRaster.Finally, we would like to thank two anonymous reviewers for their constructive comments, which improved the quality of the manuscript.

Fig. 3 .
Fig. 3. Hourly hydrograph for the Upper Ourthe at Tabreux: observed (solid blue line); HBV-96 simulation using the grid-based spatially distributed version (light grey band for the 95 % confidence bounds of 64 ensemble members); HBV-96 simulation using the operational spatially lumped version (dashed black line).Histograms of the Nash-Sutcliffe (NS) model efficiency coefficient and the root-mean-square error (RMSE) between the observed and the ensemble of simulated discharges using the grid-based version of HBV-96 are shown in the upper left corner.

Fig. 4 .Fig. 4 .Fig. 4 .Fig. 4 .
Fig. 4. Five cases of the discharge observation vector y k of increasing spatial extent.Discharge gauges contained in the observation vector are plotted in black.Channel delineation using Strahler stream ordering plotted in white pixels.

Fig. 6 .
Fig. 6.Synthetic experiment.Top: Open loop model state simulations (grey lines) including the true states (black line) at location 1 in Fig. 1b.Bottom: Forecasted (grey histograms) and analysed (dashed histograms) model states at location 1 for two dates (3 November and 2 January) and considering 5 discharge observation vectors (A, B, C, D, E).The true states are indicated by asterisks and the error bars represent the 95% confidence bounds of the observation errors (ν k in Eq. 13).

Fig. 10 .Fig. 11 .
Fig.10.Real world experiment.Forecasted (grey histograms) and analysed (dashed histograms) model states at location 1 for two dates (3 November and 2 January) and considering 3 discharge observation vectors (A, B, D).The true states are indicated by asterisks and the error bars represent the 95% confidence bounds of the observation errors (ν k in Eq. 13).

Table 1 .
Catchment response time between the catchment averaged rainfall centroid and the discharge peak.

Table 2 .
The channel width corresponding to Strahler stream order number

Table 3 .
Validation at the stream gauge locations without data assimilation.The median Nash-Sutcliffe model efficiency of the 64 ensemble members for the open loop simulations and the updated simulations being assimilated using Cases A and B. Gauge numbers correspond to Fig. 1b.

Rakovec et al.: State updating of a distributed hydrological model with EnKF
advantage of a grid-based hydrological model with gridbased routing over a lumped model with a unit hydrographtype of channel delay is that the modelled discharge is represented by spatially distributed model states, which quantify the volumes of water within the channel network.This means that we do not have to explicitly consider any time O.