A past discharges assimilation system for ensemble streamflow forecasts over France – Part 1 : Description and validation of the assimilation system

Two Ensemble Streamflow Prediction Systems (ESPSs) have been set up at M ét́eo-France. They are based on the French SIM distributed hydrometeorological model. A deterministic analysis run of SIM is used to initialize the two ESPSs. In order to obtain a better initial state, a past discharges assimilation system has been implemented into this analysis SIM run, using the Best Linear Unbiased Estimator (BLUE). Its role is to improve the model soil moisture by using streamflow observations in order to better simulate streamflow. The skills of the assimilation system were assessed for a 569-day period on six different configurations, including two different physics schemes of the model (the use of an exponential profile of hydraulic conductivity or not) and, for each one, three different ways of considering the model soil moisture in the BLUE state variables. Respect of the linearity hypothesis of the BLUE was verified by assessing of the impact of iterations of the BLUE. The configuration including the use of the exponential profile of hydraulic conductivity and the combination of the moisture of the two soil layers in the state variable showed a significant improvement of streamflow simulations. It led to a significantly better simulation than the reference one, and the lowest soil moisture corrections. These results were confirmed by the study of the impacts of the past discharge assimilation system on a set of 49 independent stations. Correspondence to: G. Thirel (guillaume.thirel@jrc.ec.europa.eu)


Introduction
Improving streamflow forecasting is a key issue for preserving human lives and material, and for monitoring water resources. Much effort has been put into coupling Land Surface Models (LSMs) with hydrological models to improve the simulation of physical processes (Miller et al., 1994;Benoit et al., 2000;Habets et al., 2008), and into increasing the spatial resolution of these models. Unfortunately, not all hydrological processes are easily predictable. In particular, hydrology is very dependent on precipitation, which is a highly stochastic phenomenon, and questions remain as to capacity to supply an adequate initial state to the hydrological model.
One attempt to address the difficulty of precipitation prediction involves the use of meteorological ensemble prediction. This kind of prediction, relying mostly on meteorological ensemble forecasts forcing an hydrological model, tends to give better scores on streamflows than deterministic predictions. Much research is underway on this subject, such as the Hydrologic Ensemble Prediction EXperiment (HEPEX) (Schaake et al., 2006, and see the website http://www.hepex.org/) which "brings together hydrological and meteorological communities from around the globe to build a research project focused on advancing probabilistic hydrologic forecast techniques".
In Europe, the European Flood Alert System (EFAS) prototype (Ramos et al., 2007) is based on the European Centre for Medium-Range Weather Forecasts (ECMWF) Ensemble Prediction System (EPS) (Chessa and Lalaurette, 2001;Buizza et al., 2007) and sends alerts to European countries.
Published by Copernicus Publications on behalf of the European Geosciences Union. G. Thirel et al.: A streamflow assimilation system for ensemble streamflow forecasts over France In France, two Ensemble Streamflow Prediction Systems (ESPSs) have been set up using the ECMWF EPS (Rousset-Regimbeau et al., 2007) and the Météo-France EPS, "Prévision d'Ensemble Action de Recherche Petite Echelle Grande Echelle" (PEARP) and have been compared using statistical scores over a long period (Thirel et al., 2008).
Data assimilation combines physical and observational information on a system in order to provide a better description of the system. The benefit of data assimilation has already been amply demonstrated in meteorology and oceanography over the past decades, where it helps to provide initial conditions for numerical prediction (Ghil and Malanotte-Rizzoli, 1991). However its use in the field of hydrology is more recent. Data assimilation in hydrological modelling can be used for three main purposes: improving soil moisture states, improving streamflow predictions, and optimizing models parameters. It can be carried out by analysing soil moisture, or/and streamflow data.
For example, Reichle et al. (2002) assessed the performance of an Extended Kalman Filter (EKF) and an Ensemble Kalman Filter (EnKF) for soil moisture analysis. Rüdiger (2005) used a variational data assimilation approach for assimilating streamflows in order to retrieve root-zone soil moisture. Recently, Zaitchik et al. (2008) used GRACEretrieved soil moisture data and improved the simulation of water storage and fluxes in the Mississippi River basin. Aubert et al. (2003) developed an EnKF assimilation system for improving streamflow prediction over a Seine river subbasin. Clark et al. (2008) used the EnKF in which states in a distributed hydrological model were updated by means of streamflow observations. They demonstrated that the standard implementation of the EnKF was inappropriate because of non-linear relationships between model states and observations and that transforming streamflow into log space before computing error covariances as well as using a variant of the EnKF not requiring perturbed observations improved filter performance.
So far, few operational applications of such assimilation systems exist. Promising work was done by Komma et al. (2008) which implemented an EnKF for an Austrian basin, for real-time flood forecasting. This system adjusts soil moisture for better real-time streamflows forecasting. Seo et al. (2009) give details of an operational variational assimilation (VAR) of streamflow, precipitation and potential evaporation data into lumped soil moisture accounting and routing models operating at a 1-h timestep.
This paper presents the work performed using assimilation to update soil moisture states of the Météo-France hydrometeorological model SAFRAN-ISBA-MODCOU (SIM), in order to improve streamflow predictions. Various soil moisture states and soil water physics are assessed in this framework. The originality and difficulty of this study lies in the fact that the data assimilation system is applied over a distributed model, for embedded station networks, and for all of France. The SIM hydrometeorological model is described in Sect. 2. Then the coupler software PALM, in which the assimilation system was implemented, and the BLUE assimilation method are described in Sect. 3. Section 4 presents the design and methodology of the assimilation system. The results obtained by the assimilation system on the SIM analysis suite with different settings are presented and discussed in Sect. 5, and a summary and a conclusion are given in Sect. 6.

The SAFRAN-ISBA-MODCOU hydrometeorological model
Both the analysis suite and the hydrological forecasts are based on the SAFRAN-ISBA-MODCOU hydrometeorological suite. This suite is composed of three independent models: SAFRAN, ISBA and MODCOU. SAFRAN (Système d'Analyse Fournissant des Renseignements Atmosphériquesà la Neige, an analysis system that provides atmospheric data to a snow model) is a near-surface meteorological analysis system (Durand et al., 1993). It combines meteorological model outputs with surface observations to produce hourly values of meteorological variables. SAFRAN provides eight parameters (10-m wind speed, 2-m relative humidity, 2-m air temperature, total cloud cover, incoming solar and atmospheric/terrestrial radiation, snowfall and rainfall) interpolated over France on the ISBA 8-km grid. Recently, Quintana  assessed the quality of SAFRAN against observations, showing that most of the parameters are well reproduced by SAFRAN.
ISBA (Interactions between Soil, Biosphere and Atmosphere, Noilhan and Planton, 1989;Noilhan and Mahfouf, 1996) is an LSM developed at Météo-France. It simulates water and energy fluxes between the soil and the atmosphere ( Fig. 1) with a simple parameterization. ISBA is used in research, numerical weather prediction and climate modelling at Météo-France. For hydrological applications (i.e. the SIM suite), the three-layer force-restore version is used (Boone et al., 1999) together with an explicit snow model (Boone and Etchevers, 2001) (Fig. 1). A subgrid runoff scheme (Habets et al., 1999a) and a subgrid drainage scheme (Habets et al., 1999b) have been implemented to tackle the issue of physical processes occuring at smaller scales than the 8-km ISBA grid. ISBA simulates the runoff through the Dunne mechanism over saturation. For soil moisture under the saturation point, the subgrid runoff is activated, its amount being lower below the field capacity, and zero below the wilting point. Drainage is produced for soil moisture above the field capacity, and residual drainage is effective below this value where no aquifer layer is present in MODCOU (see Quintana Seguí et al. (2009) for more details about the runoff and drainage processes). Recently, Quintana Seguí et al. (2009) introduced an optional exponential profile of hydraulic conductivity in the soil into ISBA, resulting in a better simulation of river discharges. This feature is intended to reduce the drainage flux to the river, which was too high at early times after heavy rains, by spreading the flux over time. MODCOU (MODèle COUplé (Coupled Model), Ledoux et al., 1989) is a distributed hydrogeological model. It simulates the spatial and temporal evolution of two aquifers (located over the Seine and Rhône basins) using a diffusivity equation. The interaction between these aquifers and the rivers is described, and the soil water is routed towards and into the rivers with a simple isochronism algorithm. Streamflows are produced with a 3-h time step, but used and validated at a 1-day time step. The ISBA drainage and runoff variables are used by MODCOU in the SIM suite (corresponding to drainage and runoff variables in italics in Fig. 1).
SIM was first validated for three large French river basins: the Rhône , the Adour-Garonne (Morel, 2003) and the Seine (Rousset et al., 2004). Then, SIM was extended and validated over the whole of France , supplying realistic water and energy budgets, streamflows, aquifer levels and snowpack simulations. Around 900 streamflow stations are simulated over France. SIM has been running operationally once a day at Météo-France since 2003 in an analysis mode. It is used for soil water reports and as a tool for the French national flood alert services, for both its streamflow and soil moisture outputs.
Based on the SIM suite, two ensemble hydrological forecast systems have been built, using the ECMWF EPS (Rousset-Regimbeau et al., 2007) and the PEARP EPS (Thirel et al., 2008). The initial soil, river and aquifer states of these two systems come from the operational analysis SIM suite described above. However, this suite is not perfect and the error in the precipitation data or in the estimation of model fluxes can lead to a bad estimation of the current state. The impact of the quality of the meteorological forecasts has already been assessed (Thirel et al., 2008), and studies on the model are ongoing, but none has been performed so far on the initial states of the model. That is why a streamflow assimilation system has been set up in SIM, in order to improve streamflow predictions. The assimilation system will rely on modifying the soil moisture of ISBA because this variable is very relevant to the river flow in the medium term. Directly modifying the amount of water in the rivers would only tackle the short term and modifying the aquifer layers would only concern the Seine and Rhône basins. The assimilation system and its impacts on the SIM suite forced by analysed data are described in the following section.

Tools used for the data assimilation system
The streamflow assimilation system was implemented in the PALM coupling software, and a linear estimation (BLUE) was used to optimize the ISBA soil moisture.

The PALM coupling software
PALM (Parallel Assimilation with a Lot of Modularity; Lagarde et al., 2001) is a dynamic parallel coupler implemented by the CERFACS (European Centre for Research and Advanced Training in Scientific Computation). PALM was written because the CERFACS was given the task of designing software that could handle the numerous methods of data assimilation needed for the oceanographic project MERCA-TOR (Brasseur et al., 2005). The specificities of PALM are a dynamic launch of the coupled components, independence of the various components which allows full modularity, and a set of standard algebra libraries (Fouilloux and Piacentini, 1999;Buis et al., 2006). Moreover, PALM is particularly well adapted to the Météo-France NEC supercomputer platform, and takes advantage of its cluster structure requiring little parallelization knowledge from the user.
All the above reasons led us to choose PALM for the implementation of the streamflow assimilation system.

The Best Linear Unbiased Estimator (BLUE) method
The BLUE method is the analysis operator used for the streamflow assimilation system. This method assumes that background and observation errors are unbiased and noncorrelated. The analysis state x a is an estimation of the true state x t , such that x a = x t + ǫ a , where ǫ a is the analysis error. The BLUE relies on minimizing T r(A) with respect to K (Bouttier and Courtier, 1999), with A, the analysis error covariance matrix, defined as follows: where R and B represent the observation and background error covariance matrices respectively, H is the Jacobian matrix of the observation operator H computed around the background state x b , and K is a linear operator (the gain matrix) to be defined. Bouttier and Courtier (1999) showed that: which highlights the fact that the data assimilation system gives a correction applied to the background state x b . y 0 is the observation vector and H(x b ) the model equivalent of the observations. The value of K minimizing T r(A) is: The H, R, y 0 and H(x b ) quantities can contain information at several time steps, depending on the size of the assimilation window. In that way, the BLUE analysis tries to find the state at the beginning of the assimilation window (x a ) that will result in the simulation closest to a set of available observations, given an a priori state (x b ) for this initial condition. The BLUE method was chosen because of the small size of the observation and state variables, which made it possible to compute the exact solution for the K matrix. It relies on the assumption that the operator H is not too non-linear over the [x b ,x a ] interval.
In the case of our application, the streamflows are assumed to be inexact, mostly because of soil moisture errors in ISBA. Thus the ISBA soil moisture is chosen to be the variable state (x b and x a variables) of the optimization process. The observations (y 0 ) used to correct soil moisture errors are streamflow observations. Consequently, B represents ISBA soil moisture error statistics and R represents streamflow observations error statistics. H(x b ) stands for streamflows computed by the SIM suite using the background soil moisture x b . H represents the model suite ISBA-MODCOU and H is its tangent linear version, computed around a reference often chosen as x b .

Streamflow assimilation methodology
The originality of the present streamflow assimilation system is that it is applied to a distributed hydrometeorological model over the whole of France. Therefore, a wide range of basins (large or small, contrasted or not) and meteorological conditions are encompassed. Moreover, single basins and embedded basin networks are assimilated simultaneously. In the following, we show how some difficulties associated with these features were overcome.

Principle of the assimilation process for SIM
The principle of the assimilation process is shown in Fig. 2 for an N-day time window, initializing ensemble forecasts lasting P + 1 days and beginning on day (D).
The background state x b is the initial ISBA soil moisture state at day (D − N). The innovation vector is the difference between the streamflow observations and the simulated streamflows (resulting from a SIM run initialized by the ISBA soil moisture at day (D − N)), from day (D − N) to day (D − 1). The observation error covariance matrix, R, is computed from statistics on the [D − N;D − 1] streamflow observations, and the background covariance error matrix B is taken to be constant over time (more details later in Sect. 4.4.1). The linear approximation of the Jacobian matrix H is computed around the background state x b with small perturbations.
Then, BLUE uses all these elements to identify the analysis soil moisture state. The (D−N) background soil moisture is corrected by this analysis state, and SIM is integrated over the [D − N;D − 1] time window. This integration provides soil moisture and river initial states for performing ensemble streamflow forecasts from day (D) to day (D + P ).
The assimilation process can be re-iterated after a delay of at least N days so as not to use the same streamflow observations several times.

Selection of gauge stations for observations
The streamflow observations come from the data collected by the "Banque Hydro" over a network of approximately 3500 river gauge stations . This French database is available online at the following website: http://www.hydro.eaufrance.fr/. The stations simulated by SIM (≈ 900) all correspond to real gauge stations present in the "Banque Hydro" database. A set of 186 relevant gauge stations (good quality of streamflow measurements and SIM results) was selected for assimilation. Thus the size of the variable state vector was 186. The discharge observations are available daily, and are dailyaveraged values.
In order to respect the river structures and to deal with dependencies between sub-basins, all the stations were sorted into river "trees" (main basins) in which the bases were the down stream station (base station) and the "branches" were its upstream sub-stations (see an example in Fig. 3). The number of stations in a tree ranged from one (independent basins) to 34 (the Loire main basin).

State variable definitions
The state variable of the assimilation system is the soil moisture of SIM. In SIM, which includes the ISBA-3L version of ISBA (Boone et al., 1999), the soil is divided into three layers: a thin surface layer, a root layer, and a deep layer (see Principle of the assimilation process for SIM for an assimilation over a N -day time window, initializing an ensemble streamflow forecasts lasting P days.

Fig. 3.
Simple representation of a basin with its three internal independent sub-basins. The subscript 1 is for the base sub-basin/station streamflow, 2 and 3 stand for the upstream sub-basins/station streamflows. y i is for the streamflow at the gauge station i, and x i represents the variable state (calculated according to the chosen method), excluding the surfaces belonging also to a sub-basin located upstream of it (i.e. x 1 only includes the non-shaded area, x 2 the horizontally shaded area, and x 3 the vertically shaded area).
of the root layer and has no impact on streamflows (this layer is mainly used to determine the surface humidity for bare soil evaporation), so the root-and deep-layers moisture are the only relevant state variables. Three definitions of the state variable were considered. If w 2 and w 3 (a different value for each ISBA grid mesh) stand for the soil water content (in m 3 /m 3 ) of the root and deep layers, respectively (see Fig. 1), and d 2 and d 3 − d 2 are their corresponding thicknesses, the elements of the first state variable are: The second method only used the root-layer water content: The last method considered each soil layer water content separately: The size of this last state variable is twice that of the previously described state variables. Each that the sum is taken over the sub-basin i, excluding meshes belonging to another upstream assimilated sub-basin.
Because the BLUE analysis provides a correction over each sub-basin i, not over each of its ISBA meshes, the final step of the assimilation process was to disaggregate the soil moisture from the analysis state space (186 values) to the ISBA grid. A simple comparison between elements of x a and x b gives a coefficient to be applied uniformly to the ISBA grid soil moisture over all meshes of the corresponding basin. The formula is given in Eq. (7) for an ISBA mesh included in the i th assimilated basin (for Eqs. 4 or 5):

Practical implementation
In the following, equations and matrices are illustrated for the simplest case of a 1-day assimilation window, and for a state variable defined with Eqs. (4) or (5) only. The equations and matrices can be easily generalized for a longer assimilation window, or for the case where soil moisture is taken into account separately for layers 2 and 3 (Eq. 6).

Specification of background and observation errors
In Eq.
(2), the background (B) and observation (R) error covariance matrices are the terms that define the modelled soil moisture and streamflow observations error statistics. So their specification is a key point of the assimilation process. Both matrices were taken to be diagonal, in order to simplify this first study. This means that the error on soil moisture for a given sub-basin was assumed not to be correlated with the error on soil moisture of any other sub-basin, and that the error on an observed discharge was not correlated with any other observed discharge. It will be demonstrated below that such an assumption does not prevent the system from being efficient. The variance of the background error was estimated by applying a known error to SAFRAN precipitation and temperature (consistent with the findings of Quintana , and examining the resulting error on ISBA soil moisture. The variance error was estimated by comparing the soil moisture of a reference SIM run with that of a SIM run forced by a perturbed SAFRAN temperature and precipitation. The two parameters were perturbed over a period of 19 months (from March 2005 to September 2006) by Gaussian white noise with rmse around 1.5 • C for temperature, and a noise with rmse around 2.4 mm day −1 for rainfall (values taken from Quintana . The variances of the background covariance error matrix were computed according to the definition of the variable state. The B elements had a mean around 10% of the square of an averaged soil water content (0.25 m 3 /m 3 ), depending on the chosen variable state.
The variance of the observation error was defined using the quantiles 1 (Q 1 ) of streamflow observations (daily flow that is exceeded 99% of the time as provided by the "Banque Hydro" database). For streamflows under this quantile, the observation variance errors were defined to be proportional to Q 2 1 (i.e. the errors on measurements were proportional to Q 1 ), and above Q 1 they were taken about (7%) 2 of the square of the streamflow observations (corresponding to measurement error proportional to 7% of the measured streamflow). This method was chosen in the following, after being compared with another method.

Jacobian of the observation operator
The observation operator H describes the link between the variable to be improved (the simulated streamflows y) and the state variable (the soil moisture x). In Eq. (2), H, called the Jacobian matrix, is the linear approximation of H and can be written (on x = x b ): Assuming the validity of the tangent-linear hypothesis, the modelled streamflow consecutive to a variation x of the initial soil moisture can be approximated by: So that, using an uncentred finite difference scheme, we have: y i is the modification of the sub-basin i streamflow resulting from a modification x j of the sub-basin j soil moisture. The computation of H consists of comparing the perturbed response of the MODCOU streamflows to a reference simulation of SIM.
However, since assimilated sub-basins are embedded in larger basins, a single perturbed run of SIM is not enough to deduce all the elements of H. In a given basin, all the subbasins have to be perturbed separately, in order to deduce the specific influence of each sub-basin on all its down stream gauge stations discharges. The detailed computation of H is given for a simple theoretical example in Appendix A.
The underlying linearity hypothesis used to derive the BLUE equation imposed the use of SIM, during the assimilation process, in domains where the model remained almost linear. To check that such an hypothesis was satisfied, a sensitivity study was performed. A range of perturbations ( x j from 0 to 10% of the initial soil water content x b ) was tested and showed that, for an applied perturbation of around 0.1%, the values of the Jacobian matrix coefficients were nearly constant with the applied perturbation. Moreover, it was shown that an opposite perturbation (−0.1%) led to a similar Hydrol. Earth Syst. Sci., 14, 1623-1637, 2010 www.hydrol-earth-syst-sci.net/14/1623/2010/ Jacobian matrix. Therefore, in the assimilation experiments presented below, all the Jacobian matrices were computed with a 0.1% perturbation applied to the soil water content.
Because of soil moisture heterogeneities in space and time, the Jacobian matrices were recalculated for each assimilation window.

Twin experiments
Several twin experiments (experiments based on synthetic observations) were performed in order to validate the assimilation system. Such experiments allow the behaviour of the assimilation system alone to be evaluated. The experiments were carried out over a 3-month period (17 February 2006 to 17 May 2006) characterized, for most of the gauge stations, by several flood events during the first half, followed by a dryer period. The analysis state included the soil water content of both root and deep zones of ISBA soil (Eq. 4).
The assimilation was performed every 5 days on a 5-day time window, with the standard physics of ISBA.
The tests consisted of adding or removing 5% to 10% of the soil water content of all the assimilated basins into the initial state of ISBA soil moisture on 17 February 2006 (4 different experiments, Table 1). These initial conditions induced severe floods or droughts, at the beginning of the simulated period, which tended to diminish over time. A reference SIM run was used to generate the synthetic observations.
The results of the 4 assimilation studies are presented in Table 1 and compared to 4 perturbed experiments without any assimilation performed. Table 1 presents the Ratio-rmse (see Appendix B for definition) of simulated SIM streamflows consecutive to the initial perturbation, during the first five assimilation periods. The first experiment (+10% initial state) shows that the assimilation was very effective, reduc-ing the Ratio-rmse to values under 0.15 even for the first assimilation time step, while this score was largely higher than 1 without assimilation for the first assimilation periods. The second experiment (+5% initial state) had the same global behaviour. However, the −5% and −10% experiments did not behave in the same way. The first assimilation time step for these two experiments seemed to be useless, with a Ratio-rmse of the same order as in the corresponding nonassimilated experiment or even higher. In fact, the assimilation "over-corrected" the error. Then, the following assimilations reduced the Ratio-rmse to values lower than 0.1, which was markedly better than the non-assimilated experiment. This "over-correction" (the first assimilation time step led, in fact, to a soil moisture 1 to 2% wetter on average than the reference state) was probably due to the non-respect of the linearity hypothesis: the H matrix was computed (for the first assimilation process) for dry values of soil moisture (largely below the field capacity value w fc ) due to the perturbed initial condition. The non-perturbed values (i.e. the "truth") of soil moisture for the initial state were around or above the field capacity value. Discharges are more dependent on soil moisture for wet soils than for too dry soils. So, since the behaviour of the physics was rather different between the background state and the analysis state, the linear hypothesis was not respected for this first assimilation. However, the system converged rapidly, and it seems important to note that the initial perturbations imposed for these twin experiments were unrealistic and, indeed, huge. For real cases, as described in the following, the increments given by the BLUE are smaller.

Assimilation of real observations
Six experiments are described in this section: the three previously described state variables were used (see Eqs. 4, 5 and 6) and, for each one, two different physics schemes were tested in ISBA (with or without the exponential profile of hydraulic conductivity in the soil). These experiments are summarized in Table 2 with a reference name that is used in all that follows. For the IS 1 , IS 3 and IS 5 experiments, the balance between R and B matrices was chosen by testing a range of possible balances on the shorter 3-month period previously used. Then this was extended to, respectively, IS 2 , IS 4 and IS 6 . The goal of this comparison was to find the best possible set of initial states for the ensemble streamflow forecasts based on SIM. That is why the chosen study period (from 10 March 2005 to 30 September 2006) corresponded to the period for which the comparison of the two ensemble streamflow forecasts chains of Météo-France had been performed by Thirel et al. (2008). However, the assimilation was started on 2 January 2005 in order to allow a 2-month spin up of the system.
The assimilation of observations was done daily in order the system to react to bad simulations fast. Therefore, the assimilation window was 1-day long. In order to limit the increments not respecting the validity of the linear hypothesis of the BLUE, the increments were limited to a ±10% range. The validity of this hypothesis will be assessed in Sect. 5.4. Moreover, for dry soils (soil moisture lower than 1.1w fc ), negative increments were limited to −2% since during observed low flows, if the SIM streamflow was overestimated during several consecutive days, the BLUE tended to dry the soil moisture down to very weak and unrealistic values. This tendency did not actually improve streamflow simulations, and resulted in a severe underestimation of the first few following flood events. This modelling problem can be explained by bad simulation of an aquifer, high anthropization of the basin or difficulties in measurements, which cannot be resolved by adjusting the ISBA soil moisture. Figure 4 presents the accumulated distribution of efficiency for the 186 assimilated stations, for the 6 experiments plus the two reference simulations of SIM. It shows that the best simulations are IS 1 , IS 2 and IS 5 , and that the improvement in the Nash criterion is significant. For each value of efficiency, the three solid lines are largely higher than the reference solid line. This shows that the experiments with the standard physics of the model (no exponential profile of the hydraulic conductivity) are better than the no-assimilation reference, especially for the IS 1 and IS 5 experiments. The dashed lines (experiments with the exponential profile of the hydraulic conductivity) are closer to the reference dashed line, because the physics improves the streamflow simulation, but remain above it.
The scores (see Appendix B for definitions) presented in Table 3 are averaged for a selection of 148 stations all over France, out of the 186 available (as shown in Fig. 5, left). Stations for which the data assimilation system did not improve the streamflow simulation were excluded here. These stations are located on down-stream parts of the Seine or the Rhône aquifer layers (which are explicitly simulated by SIM). For these basins, the streamflow is mostly affected by the aquifer level rather than by the precipitation or the Table 2. Definition of the six different experiments assessed for the assimilation of real observations. IS stands for "Initial States", because these states will be used as initial states for the Ensemble Streamflow forecasts in the future.

Experiment
State variable Exponential profile Eq. (4): w 2 + w 3 Yes IS 3 Eq. (5) soil water content. Moreover, for some stations, the missing observations were too important so these stations were also excluded. Table 3 shows that, for every experiment, the discharges simulations were improved (better Nash criterion and Ratiormse than reference). The assimilation system for experiments without the exponential profile of hydraulic conductivity (IS 1 , IS 3 , IS 5 ) allowed a greater improvement of the streamflow simulation (when compared to the reference) than the IS 2 , IS 4 , IS 6 experiments. Since the reference simulation was less accurate here, more correction could be made Hydrol. Earth Syst. Sci., 14, 1623-1637, 2010 www.hydrol-earth-syst-sci.net/14/1623/2010/ . Nash, Q, Ratio-bias and Ratio-rmse scores are presented for each assimilation experiment (straight) and its corresponding reference simulation (italics). "Mean incr. (abs.)" stands for the mean of the absolute values of the increments imposed by the BLUE analysis, and "Incr. σ " stands for the spread of these increments. As the BLUE gives separate increments for layers 2 and 3 of ISBA for experiments IS 5 and IS 6 , the two last columns are split into two sub-columns with mean and spread of increments for layer 2 (left) and for layer 3 (right  by the assimilation. Moreover, it can be seen that the Q ratio was acceptable for all experiments and that assimilation always reduced the model bias. The Ratio-rmse, the mean of absolute values of increments, and the spread of the increments were lower for the experiment with the exponential profile of hydraulic conductivity (with the same state variable). This means that the assimilation system performed better for this experiment (i.e. the Ratio-rmse was reduced), and that smaller changes were imposed by the assimilation system, that is to say the water budget of SIM was less modified. This indicates that the exponential profile of hydraulic conductivity should be chosen, rather than the standard physics of the ISBA scheme, in order to limit the modification of the ISBA prognostic variables by the assimilation. Nevertheless, experiment IS 5 had a better Nash criterion than experiment IS 6 , for which the only difference was the use of an exponential profile of hydraulic conductivity. This poor performance of the IS 6 experiment can be explained by the fact that the assimilation system acted, in an independent way, on the third soil moisture layer. So it acted directly on the drainage flux to the river, which was a phenomenon having a much longer time-scale than the 1day assimilation window, with the exponential profile of hydraulic conductivity version of ISBA.
It clearly appears that the state variable defined by Eq. (4) (average of the root-and deep-layer soil moistures) gave better results than when defined by Eq. (5) (root layer only). The Nash criterion was the best for these experiments, and the Ratio-rmse, mean increments and spread of increments were lower. The superiority of the Eq. (4) definition over the Table 4. Statistical scores for SIM discharges and the BLUE assimilation system for a set of 49 independent stations. Scores are averaged over a 19-month period (10 March 2005-30 September 2006. Nash, Q, Ratio-bias and Ratio-rmse scores are presented for each assimilation experiment (straight) and its corresponding reference simulation (italics). Assi.

Ref.
Assi. Ref. Eq. (5) one can be explained by the fact that changing only the root-layer soil moisture (i.e. Eq. (5) definition) controls the streamflow simulation driven by the runoff alone, not by the drainage. Moreover, a drainage flux of soil water content goes from the root layer to the deep layer, modifying the streamflow simulation, and there is no possibility for the assimilation system to directly act on the deep layer if too much water is added there. With a comparable Nash criterion of about 0.8, the three best experiments were IS 1 (Nash = 0.81), IS 2 (Nash = 0.80), and IS 5 (Nash = 0.83). Despite a better Nash criterion, the IS 5 experiment showed higher Ratio-bias, Ratio-rmse and increment values than IS 2 , revealing less efficient behaviour. Moreover, because the physics of IS 2 was improved with respect to the physics of IS 5 , it can be assumed that the assimilation system contribution should be maintained over longer periods. It is also important to notice that the CPU time cost of IS 5 is twice that of IS 1 or IS 2 . The same conclusions (except for the computational cost) can be drawn between IS 1 and IS 2 , illustrating that IS 2 should be preferred to IS 1 .
The best configuration for assimilating the streamflows in SIM appears to be IS 2 . Moreover, the mean of the absolute values of the increments (0.23%) for this experiment is comparable to the value, 0.1%, of the perturbation used to compute the Jacobian matrix for the BLUE. This reveals that the assimilation system should respect the linearity hypothesis of the BLUE most of the time. This will be assessed below.
The behaviour of the IS 2 assimilation process during a shorter 200-day period is examined in Fig. 6 for streamflows (Fig. 6a), layer 2 and 3 soil moistures (averaged over the subbasin considered, Fig. 6b), and BLUE increments (Fig. 6c), for the specific case of the River Doubs at Besançon. This period was characterized by a wet period with several flood events, followed by a dryer period. The assimilation was remarkably efficient for the flood events. For the largest flood, the assimilated streamflow was very close to the observation (560 m 3 s −1 ) whereas the no-assimilation simulation did not produce values beyond 320 m 3 s −1 . The assimilation system seemed quite sensitive and reactive during this period, with mean absolute values of increments being regularly larger than 1% of the soil water content. The soil moisture was wetter with assimilation than without. For the largest flood, the increment was very large (+8%), allowing SIM to immediately simulate a higher streamflow than it would have otherwise.
The assimilation performed fewer corrections during the following dryer period. In fact, hardly any increments were produced during this period and the streamflows were lower than the observations, with few improvements when compared to the reference SIM simulation. This poor performance of the assimilation was caused by the fact that the soil was rather dry. Because the soil moisture was lower than the field capacity, streamflows had a lower dependence on the soil moisture, the model and rainfall forcing characteristics being more important for this case.

Validation over independent stations
All the scores previously presented were computed over stations whose discharge observations had been used by the assimilation system. Such a validation is not sufficient to prove that the assimilation works correctly (Talagrand, 2003). For this reason, a selection of 49 independent station was used for another validation. These stations, not used by the assimilation system, were located upstream or downstream of the stations used in the assimilation system (see Fig. 5, right). They thus benefited from the better soil moisture of the subbasins concerned, and from the improved river water content.
Some scores concerning this independent validation are presented in Table 4. Although the Nash, Q, Ratio-bias and Ratio-rmse are better for these 49 independent stations than for the 148 stations studied previously when we look at the experiments without assimilation, the scores with assimilation do not show the same behaviour. The 148 stations whose observed discharges were used by the assimilation system have better scores, which is logical. However, Hydrol. Earth Syst. Sci., 14, 1623-1637, 2010 www.hydrol-earth-syst-sci.net/14/1623/2010/ it is very interesting that the scores for the 49 independent stations were improved by the assimilation of streamflow observations (except for the IS 4 experiment). This overall improvement of the scores of the experiments using the assimilation system shows that this assimilation system is actually effective for improving the SIM river flow simulation.
Here again, the experiments IS 3 and IS 4 showed the lowest scores. This was due to the fact that these experiments were already the least efficient for the 148 used stations, and also to the fact that the increments of these two experiments were the highest of the 6 experiments. Thus, with high increments, fluxes were more modified, and the adjustments of the soil moisture did not necessarily lead to discharges simulations fitting the non-used observed discharges. The scores of the four other experiments were very close, with Nash criteria between 0.70 and 0.74, a Q ratio of 1, a Ratio-bias close to zero, and a Ratio-rmse between 0.61 and 0.64. However, for each of these scores, the IS 2 experiment was the best, confirming the conclusions of the study of the 148 dependent stations' scores. Although the initial states are of a quite similar quality, the use of a better physics scheme, combined with small modifications of the soil moisture, are key ingredients for an improvement of discharges on independent stations. That is why this configuration (IS 2 , which uses the improved physics and an assimilation method modifying the soil moisture of the deep and root layers jointly) will be chosen to initialize the ensemble streamflow forecasts (Thirel et al., 2010).

Impact of iterations of the BLUE
Using the BLUE to implement a data assimilation system requires using the model in an area of linear behaviour. The linearity is necessary for the BLUE to provide the exact solution to the analysis equation. Of course, with many models, and particularly with hydrological models, linearity is not always respected. It was therefore necessary, in order to verify the relevance of using the BLUE, to investigate this point.
An external loop was implemented on the BLUE, i.e. the BLUE was iterated, until convergence, around the analysis state. In practice, this means that a first iteration of the BLUE was performed (as previously) but then another iteration was performed by computing a new Jacobian matrix around the analysis state given by the first iteration. For this second iteration, the background state modified by this new iteration of the BLUE was once again the one used for the first iteration. Such iterations were performed until the decrease of the Ratio-rmse was lower than 0.1 for all the assimilated stations, and a maximum of 10 iterations was fixed. It is obvious that such a process needs at least twice as much CPU time (when the convergence is immediate) and possibly 10 times as much CPU time (the maximum number of iterations). Thus, testing of the iterations for this study was limited to the IS 2 experiment.
The results of the iterations of the BLUE on the scores of the IS 2 experiment are presented in Table 5 together with the results of SIM without assimilation but with the improved physics, and the results of the original IS 2 experiment. This table shows a very slight improvement of the Nash criterion (from 0.80 to 0.81) and of the Ratio-rmse (from 0.53 to 0.51) when the iterations of the BLUE are used instead of the original BLUE. However, the Q ratio and the Ratio-bias are slightly deteriorated by the use of the iterations of the BLUE. Finally, the mean of the absolute values of the increments is not modified. All these scores indicate that the contribution of the iterations is very weak, so, because of its prohibitive CPU time cost, its use is far from worthwhile. Moreover, this experiment showed that the linearity of SIM was quite well respected by the past discharge assimilation system, justifying the use of the BLUE for this particular application.

Conclusions
The implementation of a streamflow assimilation system in the SIM distributed hydrometeorological model has been described in this paper. The performance of this system was assessed for four twin experiments, which showed that the system was efficient in the case of assimilation of synthetic observations. Then the system was assessed for six different configurations over a 569-day period with a set of statistical scores for the assimilation of real observations on 148 dependent stations discharges, and also on 49 independent station discharges. Finally, the respect of the linearity of the model was assessed by studying the contribution of an external loop.
The BLUE theory was used for the data assimilation system in the modular coupler PALM. The assimilation system was designed to correct the model soil moisture in order to bring the simulated streamflows closer to their true state. Streamflow observations from a total of 186 gauge stations were used. Three choices for the state variable were assessed: considering both layer 2 and layer 3 soil moisture in a single state variable, only considering the layer 2 soil moisture, or considering the soil moisture of the two layers separately. Moreover, the impact of using an improvement in the physics (the use of the exponential profile of hydraulic conductivity) was assessed for each of these choices.
This assimilation system showed results that were very encouraging for the application of such a method to the SIM distributed hydrometeorological model. An overall improvement of the Nash criterion, and a reduction of Ratio-bias and Ratio-rmse were observed for each configuration considered, for a selection of 148 assimilated stations. It is important to note that the use of the exponential profile of hydraulic conductivity did not necessarily improve the Nash criterion but, as it reduced the Ratio-rmse, as well as the size of the increments produced by the assimilation system, this physics should be used. For an equivalent streamflow simulation, the less the soil moisture is modified, the better the system performs. Moreover, the benefit of the improved initial state of the assimilation would last longer with this better physics.
The best assimilation (combining high scores and low modification of the model soil moisture) was found for a state variable that was the mean of soil moistures for the two layers with the exponential profile of hydraulic conductivity (IS 2 experiment). When only the root layer was used, increments were seen to be too large and poor performance was noted. Adjusting the two layers separately could give the best Nash criterion (IS 5 experiment) but the overall performance of the assimilation system (bias, rmse, increments) was lower than for the IS 2 experiment. The IS 2 configuration should thus be chosen, in order to combine good performance of the simulation, small soil moisture corrections, and the respect of the Q ratio.  . "No assimilation" was the experiment without data assimilation but using the improved physics, "IS 2 " used the data assimilation system, and "IS 2 with iterations" was the data assimilation experiment using the BLUE iterations until convergence. The three experiments used the improved physics. These conclusions were reinforced by the study of the scores for 49 independent stations. These scores showed the best behaviour of the IS 2 experiment for the Nash criterion, the Q ratio, the Ratio-bias and the Ratio-rmse. This was due to the better physics used in this configuration and to the low increments imposed by the BLUE.
Finally, it was shown that the hypothesis of the linearity of H was quite well respected when using the BLUE, as the use of an external loop increased the performance of the system only a little. As the CPU time cost of an external loop is very high, it will not be selected. The improved physics greatly reduced the intensity of the increments, as, it did the respect of the hypothesis of the linearity of H, showing the interest of using it.
Because of the lack of efficiency of adjusting the soil moisture for basins where the aquifer layers simulated in the SIM model ruled the streamflows, an assimilation of aquifer levels should be combined with the present assimilation system. Moreover, a calibration of the balance between R and B for each assimilated basin could improve the performance of the SIM model. The modularity of the PALM coupler, and the structure of the algorithm as implemented in PALM (perturbed runs of SIM) could also allow the estimation of K in the BLUE to be easily replaced using Ensemble Kalman Filter, an approach that could handle non-linear effects of the model better (even though these non-linear aspects seemed negligible). Finally, the assimilated basins could be reorganized in order to save CPU time. It could be interesting to subdivide the biggest basins into a smaller number of assimilated sub-basins, thereby reducing the number of perturbed runs of SIM needed to compute the Jacobian matrix.
This study is specific for the SIM model and the France area. The assimilation system has to be adapted to the model used and its variables, and to the area-specific hydrological variables. For example, assimilating rainfall radar observations would surely prove a better efficiency for arid areas, but for snowfall dominated areas, the modification or assimilation of the snowpack would be more relevant.
This study has demonstrated the potential of using a past discharge assimilation system in order to improve the SIM streamflow simulations and then to provide good quality initial states for ensemble prediction systems. The impact of the IS 2 initial states on the scores of two Ensemble Streamflow Prediction Systems of Météo-France (PEARPand ECMWF-based SIM hydrological forecasts) will be examined in a forthcoming study (Thirel et al., 2010).

Filling of the Jacobian matrix for a simple main basin
A simple, theoretical example is given here to explain how the H matrix was computed. Figure 3 represents a schematic river network, where three stations are assimilated. Two upstream sub-basins (subscripts 2 and 3), independent of each other, flow into a larger down-stream sub-basin (subscript 1). x stands for the soil moisture as the state variable, and y stands for the streamflow as the observation.
With these three assimilated stations, the Jacobian matrix is a 3 × 3 matrix (for a "tree"): Three perturbed runs of SIM (in addition to a reference run) had to be processed to estimate H completely. For the first one, only basin 1 soil moisture was perturbed, with no change to the soil moisture of sub-basins 2 and 3 . With this run, the response of station 1 streamflow (y 1 ) to a perturbation of the soil moisture of sub-basin 1 (x 1 ) was deduced, and so H 1,1 was known (the state variable element corresponding to sub-basin 1 excludes the soil moisture of sub-basins 2 and 3, so its change has no effect on streamflows 2 and 3 (i.e. H 2,1 and H 3,1 were null)). Then, when x 2 was perturbed, the response of y 2 (H 2,2 ) and y 1 (H 1,2 ) were simultaneously deduced. And finally, a last run, with a perturbation on x 3 , allowed the two last coefficients of H (H 1,3 and H 3,3 ) to be computed. The full Jacobian matrix for this illustrative example is displayed in Eq. (A1).
For the assimilation of the 186 stations, 34 perturbed runs of SIM were needed (plus 1 non perturbed run), because this is the number of perturbations to be imposed to compute the Jacobian elements of the largest basin (the biggest basin, the Loire basin, has 34 stations). Other independent basins can be treated during the determination of the part of the Jacobian matrix dedicated to the Loire basin.

Statistical scores
Four hydrological statistical scores were used for this study: the Ratio-root mean square error (Ratio-rmse), the Ratiobias, the discharge ratio Q, and the Nash criterion (or efficiency).
The Ratio-root mean square error (Ratio-rmse) is defined as: with T the total number of days of the period studied, Q t obs the observed discharges (for day t), Q obs the mean of the observed discharges during the study period, and Q t sim the simulated discharges (for day t).
The Ratio-bias is computed as: The discharge ratio (Q) is the ratio between the mean simulated discharges and the mean observed discharges: And finally, the Nash criterion (E, efficiency) is: The Nash criterion is a score of between −∞ and 1, and it is often considered that a score higher than 0.7 characterizes a very good simulation of the discharges. Furthermore, this score, by its definition, is very sensitive to flood errors.