Review of the Kalman-type hydrological data assimilation

ABSTRACT

There is great potential in Data Assimilation (DA) for the purposes of uncertainty identification, reduction and real-time correction of hydrological models. This paper reviews the latest developments in Kalman filters (KFs), particularly the Extended KF (EKF) and the Ensemble KF (EnKF) in hydrological DA. The hydrological DA targets, methodologies and their applicability are examined. The recent applications of the EKF and EnKF in hydrological DA are summarized and assessed critically. Furthermore, this review highlights the existing challenges in the implementation of the EKF and EnKF, especially error determination and joint parameter estimation. A detailed review of these issues would benefit not only the Kalman-type DA but also provide an important reference to other hydrological DA types.

Editor D. Koutsoyiannis; Associate editor F. Pappenberger

KEYWORDS:

• data assimilation
• extended Kalman filter
• ensemble Kalman filter
• model error
• parameter estimation

1 Introduction

In the past few decades, hydrological models have significantly benefited from the improvement of computational capacity and the availability of multi-source measurement data to simulate and forecast hydrological processes. However, due to the deepening involvement of uncertainties stemming from initial conditions, inputs and outputs, operational hydrological models will remain inherently imprecise.

Data assimilation (DA) is a procedure developed to optimally merge information from model simulations and independent observations with appropriate modelling (Liu et al. 2012b). The DA could provide optimized initial conditions, updated parameters and even improved structures for the dynamic model. Originally used in areas such as atmospheric and oceanographic science (Derber and Rosati 1989, Dee 1995), DA has drawn the attention of hydrologists more and more for its convincing performance in the real-time correction of hydrological models (Robinson and Lermusiaux 2000).

Walker and Houser (2005) and Reichle (2008) summarized the basic hydrological DA methods; Moradkhani (2008) reviewed remote-sensing measurement techniques and their applications in DA; Liu and Gupta (2007) discussed the role of DA in addressing the uncertainties in hydrological models and discussed the application of the hydrological DA in operational scenarios (Liu et al. 2012b). Montzka et al. (2012) examined the joint assimilation of observational data precedents from different spatial scales and different data types. Reviews that focus on a specific method are relatively rare compared to reviews with a wide scope (e.g. Van Leeuwen 2009).

The Kalman filter (KF) is a classic sequential method that has been widely used in hydrological DA for more than two decades (Evensen 1994b, McLaughlin 1995). Compared to other methods, such as particle filters (PFs) and variational methods, the KF is easier to implement and could produce comparable or even better results with a lower computation demand (Weerts and El Serafy 2006, Abaza et al. 2014b). Furthermore, the KF is very flexible for coupling with hydrological models and has more derivative variants than any other methods. Among these variants, the Extended KF (EKF) and the Ensemble KF (EnKF) were both developed to extend the application of the KF to nonlinear systems; also, they are the two major descendants of the linear KF. The EKF applies a straightforward Taylor extension scheme to linearize the nonlinear system. As logical as it is, the EKF has endured a “notorious” reputation for being unstable when applied to complex nonlinear hydrological models. The EnKF avoids direct linearization by statistically analysing the ensemble members. Although it increases the computational cost, the EnKF is one of the most widely used hydrological DA methods. Despite some comparison case studies (Reichle et al. 2002b, El Serafy and Mynett 2004, Dumedah and Coulibaly 2012), there is no critical review that specifically focuses on the KF. The objective of this review is to fill this gap by assessing the latest developments and analysing the challenges of Kalman-type hydrological DA, especially the EKF and EnKF. Nevertheless, this review is not intended to judge which method is better, but rather regards them as different solutions to the same problem. This is not only because they both branched from the linear KF, but also because they face some similar issues when applied in hydrological DA, of which some are also faced by methods other than KFs.

In the rest of the paper, we introduce hydrological DA, including the DA targets and methods in Section 2. The theory of the KF (including EKF and EnKF) and its state of the art applications in hydrological DA are discussed in Section 3. Section 4 discusses issues regarding the implementation of the KF in hydrological models. The summary and conclusions are given in Section 5.

2 Hydrological data assimilation

2.1 Hydrological DA targets

Data assimilation was first used in the 1950s for numerical weather forecast modelling; however, hydrologists did not pay much attention to it until the 1990s (Evensen 1994b, McLaughlin 1995). Proper use of DA may help in handling uncertainties from model inputs, the initialization and propagation of states, the model structures and even the model parameters (Vrugt et al. 2006, Liu and Gupta 2007, He et al. 2012). Meanwhile, global DA may improve regional field estimation by achieving more accurate external boundary condition estimations (Robinson and Lermusiaux 2000), and local DA could also lead to improved global estimations (Clark et al. 2008).

The development of remote sensing has promoted the application of DA in hydrological models. The remotely sensed hydrological data that have had, or currently have, the potential to be applied in hydrological models include (Houser et al. 2005, Walker and Houser 2005, Xu et al. 2014):

1. overland parameters (e.g. topography, land cover, albedo);

2. forcing inputs (e.g. precipitation, air humidity and temperature);

3. states (e.g. soil moisture, snow cover);

4. fluxes (e.g. carbon flux).

The overland parameters are usually regarded as “static” in the model, even though this is not true in long-term forecasts. Forcing inputs have the potential to replace traditional ground-based observations, with the development of remote-sensing instruments and more accurate retrieval algorithms. Neither overland parameters nor forcing inputs are the main interest of DA for hydrological modellers at this stage. The states and fluxes provide validation to the intermediate processes of the models. This was not available earlier, and hence has drawn extra attention.

Snow cover has a high albedo and thermal properties, as well as medium-term water storage capacity; therefore, the assimilation of snow observations could improve hydrological prediction (Walker et al. 2003). Andreadis and Lettenmaier (2006) used the EnKF to assimilate remotely sensed snow observation moderate resolution imaging spectroradiometer snow cover extent (MODIS SCE) data into the VIC hydrological model to update snow-water equivalent (SWE) estimates. From these estimates, a simple snow depletion curve scheme from SNOTEL station SWE data and MODIS imagery (Moradkhani 2008) was used to form the observation operator. Clark et al. (2006) assimilated the observations of a snow-covered area (SCA) to update the hydrological model, and they found that the assimilation of SCA information results in minor improvements in the accuracy of streamflow simulations near the end of the snowmelt season. Both of them found that the snow-cover update works better during the snowmelt season than the snow-accumulating season. Further research with regard to snow-cover assimilation can be found in the literature (Pan et al. 2003, Sheffield et al. 2003, Rodell and Houser 2004, Kumar et al. 2008, Parajka and Blöschl 2008).

Galantowicz et al. (1999) demonstrated KF retrieval of the soil moisture profile and temperature from L-band radio brightness observations, while Crow and Wood (2003) applied the EnKF to assimilate airborne measurements of surface brightness temperature into a TOPMODEL-based Land–Atmosphere Transfer Scheme (TOPLATS). Due to the short memory the land-surface skin temperature holds, it is suggested that it could be combined with other state variables with longer memories, such as deeper soil temperature or moisture, to obtain longer DA effectiveness (Walker et al. 2003). The surface temperature products are also sensitive to terrain, vegetation and cloud contamination; hence, multiple sets of products and better retrieval algorithms are required for continuous DA operations (Huang et al. 2008).

Soil moisture is the key variable to control the runoff generation process. It is logical to assimilate soil moisture directly, because it is a continuous storage indication as well as a natural bridge to transfer the hydrological model into state-space equations. The correct estimation of antecedent soil moisture content is critical to streamflow simulation when the soil is neither too dry nor over-saturated (Reichle et al. 2002a). Soil moisture assimilation has been proven beneficial to streamflow prediction, especially that for low flow. (Chen et al. 2014, Wanders et al. 2014b). The application of soil moisture requires the ground-based measurements to be dense enough (Chen et al. 2011), otherwise the spatial interpolation error will have to be considered. Satellite remote sensing can provide an economic and more spatial-temporal continuous observation of soil moisture.

Parrens et al. (2014) assimilated in situ soil moisture observations into a soil model at local scale. Reichle et al. (2002a) assimilated L-band microwave brightness temperature into a land surface model to estimate the near-surface soil moisture. Crow and Ryu (2009) used remotely sensed soil moisture retrievals to correct both the soil moisture state and satellite rainfall products. Alvarez-Garreton et al. (2014) assimilated surface soil moisture and soil wetness index derived from passive microwave Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E). Brocca et al. (2012) assimilated surface and root-zone soil moisture products derived from the active microwave advanced scatterometer (ASCAT) into a rainfall–runoff (RR) model using the EnKF. Recently, soil moisture assimilation has also been used to assist the parameter identification of hydrological models (Tran et al. 2014, Wanders et al. 2014a).

Despite the existence of some root-zone soil moisture by-products from the surface soil moisture observations (Wagner et al. 1999, Das et al. 2008), one of the obvious drawbacks of remote sensing of soil moisture is that only surface or near-surface soil moisture is available (Moradkhani 2008, Han et al. 2012). Preliminary experiments showed that surface soil moisture assimilation has a minimal effect on the simulation of the deep-layer soil moisture (Chen et al. 2011), although it is the latter that has a more significant impact on runoff simulations (Brocca et al. 2012, Han et al. 2012). Houser et al. (1998) argued that the remotely sensed soil moisture must be supplemented by in situ surface and root-zone observations to specify error correlation, calibrate parameters, and validate the model-calculated fields. Draper et al. (2011) found that it might be more effective to address the cause of model bias instead of relying on soil moisture assimilation to correct it. Some research combined soil moisture assimilation with the correction of other RR model forcings, such as precipitation, to improve streamflow predictions (Chen et al. 2014, Massari et al. 2014). The application of remote sensing of soil moisture may also involve issues such as rescaling (Kaheil et al. 2008, Sahoo et al. 2013), error evaluation (Alvarez-Garreton et al. 2014, Doubková et al. 2012), and radiative transfer modelling (Verhoef and Bach 2003, Reichle 2008), among others. A more detailed explanation of the limitations of soil moisture assimilation can be found in the literature (Vereecken et al. 2015).

Streamflow is the most commonly used and sometimes the only available prognostic observation variable (Clark et al. 2008, Randrianasolo et al. 2014, Samuel et al. 2014, Trudel et al. 2014, Abaza et al. 2014a). Great efforts have been made to improve streamflow forecasting using output assimilation/error assimilation over the past two decades (Broersen 2007, Anctil et al. 2003, Yu and Chen 2005, Sene 2008). The output/error assimilation methods treat a streamflow forecast as a pure model output and update it by adding errors calculated with another independent procedure or model. Such procedures/models could either be nonlinear, such as artificial neural networks (ANNs; Anctil et al. 2003), or linear, such as autoregressive–moving-average models (Broersen 2007, Chen et al. 2015). Output/error assimilation of streamflow is relatively simple to implement, as there is no feedback to the original RR model from the manipulation of model outputs.

Besides being the direct assimilation goal (Liu et al. 2012a), streamflow is also the prevalent observation to assimilate other state variables and parameters (Coustau et al. 2013). In hydrological DA, streamflow is often treated as a diagnostic variable, and is therefore not updated directly (Clark et al. 2008). In the case of a nonlinear measurement operator, it is probably the only choice to augment the state vector with the streamflow (Evensen 2003). By doing so, the nonlinear measurement operator should be reduced to a linear matrix, otherwise one would have to linearize the measurement operator. Pauwels and De Lannoy (2009) attempted to linearize the nonlinear discharge–watershed storage relationship (observation operator) in a “brutal” way. Linearization was undertaken within both a simple time series model and a conceptual Hydrologiska Byråns Vattenbalansavdelning (HBV) model; the results showed that direct linearization should be bypassed to obtain a better assimilation result.

Streamflow assimilation is different from other state variables because it involves the issue of routing. When the runoff is assimilated into a hydrological model at the current time step, not only does the current state of the watershed at a given location (e.g. outlet) need to be updated and propagated forward, but also the state at a number of different locations at multiple previous time steps (Pauwels and De Lannoy 2006). This is the case especially concerning large and distributed hydrological models. Many authors actually choose not to go into the details of the complicated channel network to control the complexity of the problem (Weerts and El Serafy 2006, Clark et al. 2008).

Despite the unique challenges faced by these targets, they share some common issues in their implementation, most notably in the quantification of observation errors and the determination of observation operators that connect them with the model output (Andreadis and Lettenmaier 2006). It is also worth pointing out that the assimilation of one single variable does not necessarily improve the estimation of the other variables. Trudel et al. (2014) reported that the assimilation of the streamflow at the outlet of the watershed would distort the estimate of soil moisture. A combination of multiple observation variables seems to be superior to single variable assimilation (Xie and Zhang 2010, Trudel et al. 2014).

The last decade saw the development of various global and regional land assimilation systems (Cosgrove et al. 2003, Mitchell et al. 2004, Rodell et al. 2004, Kumar et al. 2008). These systems enable the incorporation of multi-source observations, multi-models and multi-assimilation schemes in creating an optimal land surface states output. Due to possible water and energy balance issues, a major concern in land assimilation systems is whether they should be coupled to atmospheric models (Walker et al. 2003). Betts et al. (2003) compared the water budget of European Center for Medium range Weather Forecasting (ECMWF) 40-year reanalysis and NASA (National Aeronautics and Space Administration) DAO (Data Assimilation Office) fvGCM (finite-volume general circulation model) with the hydrological balance of the VIC model and the radiative fluxes with basin averages derived from International Satellite Cloud Climatology Project (ISCCP). They found that the runoff from both atmospheric models was significantly underestimated compared to the VIC runoff simulation, which is consistent with the observed streamflow. Large bias is also observed for radiation fluxes and surface temperature between ECMWF 40, fvGCM and the ISCCP data. Pan and Wood (2006) suggested a constrained EnKF to maintain the benefit of DA without violating the water balance principle. Boulet et al. (2000) and Bøgh et al. (2004) developed a simple water and energy balance model that allows the direct application of the remote-sensing data.

2.2 DA methods

The DA methods can be divided into different categories based on different standards (Rakovec et al. 2015). According to the dimension they focus on, these approaches can be classified as objective analysis methods and time-dependent methods (Wang and Kou 2009).

2.2.1 Objective analysis

Objective analysis aims to minimize the error between the observation field and the background field (usually the output of the numerical models) by “fusing” the new observations into the background with spatial dimensions. Typical objective analysis methods include successive correction (Cressman 1959, Barnes 1964), optimal correction (Gandin 1965, Lorenc 1981), statistical bias correction (Piani et al. 2010), Newtonian nudging (Houser et al. 1998, Paniconi et al. 2003), variational methods (Reichle et al. 2001, Seo et al. 2003, 2009, Navon 2009), etc. Detailed descriptions about the theory and development of the above methods can be found in the literature (Navon 2009).

Since objective analysis methods mostly take “snapshots” of the background filed at a given time (McLaughlin 1995), the temporal dimension is usually not well incorporated with the spatial dimensions. Even though some of them do consider the temporal dimension, such as 4-dimensional variational (4D VAR) and Newtonian nudging, time-variant objective methods are most simply viewed as a dynamic extension of the time-invariant version of the objective methods (McLaughlin 2002).

Among various objective analysis methods, the successive correction and nudging methods fail to consider the errors in the observations, while 3-dimensional variational (3D VAR) and 4D VAR ignore the uncertainties in the models. Objective analysis methods usually involve huge computations (Clark et al. 2008). Houser et al. (2005) suggested that the adjoints should be calculated as the model is developed, but this is by no means an easy task for distributed hydrological models.

2.2.2 Time-dependent methods

Time-dependent methods use a probabilistic framework and estimate the system state sequentially by propagating information forward in time (Bertino et al. 2003). The strength of the time-dependent methods is in time series analysis. The time-dependent methods work on a fixed but moving time window, and only the most recent observations that fall into this window are incorporated into the final estimation results. Typical time-dependent methods that are frequently used in hydrological DA include the linear KF (Kalman 1960) and its variants, such as the Extended KF (EKF) (Puente and Bras 1987), the Ensemble KF (EnKF) (Evensen 1994b), the Unscented KF (UKF) (Wan and Van Der Merwe 2000), the PF (Pham 2001, Moradkhani et al. 2005a, Weerts and El Serafy 2006), and the H-infinity filter (Moradkhani et al. 2005a, Wang and Cai 2008, Lü et al. 2010), among others. There are also some alternative approaches to solving specific issues in hydrological DA, most notably in the application of a genetic algorithm (GA) in the estimate of pixel-based soil hydraulic parameters for hydroclimatic modelling (Ines and Mohanty 2008) and evolutionary-based assimilation in streamflow simulations in ungauged watersheds (Dumedah and Coulibaly 2012). The spatial variation of the relationship between the background field and observation field of the variable in question is of less concern in time-dependent methods. Meanwhile, time-dependent methods are capable of handling more uncertainties (Moradkhani 2008) and are less complex to implement compared to some objective analysis methods that require an inverse or a joint model (Bertino et al. 2003).

Time-dependent methods can be well explained with Wiener–Kolmogorov estimation theory (Kolmogorov et al. 1941, Wiener 1949). Suppose is the original signal, and is the estimated signal; the estimation error is defined as: where α is the delay of the estimation. In other words, the error is the difference between the estimated signal and the true signal shifted by α. Depending on the value of α, the estimation problem can be described as in Figure 1:

1. if α > 0, the estimation is a prediction problem (error is reduced when is similar to a later value of s(t));

2. if α = 0, the estimation is a filtering problem (error is reduced when is similar to s(t));

3. if α < 0, the estimation is a smoothing problem (error is reduced when is similar to an earlier value of s(t)).

Figure 1. Types of estimation problem: (i) prediction, (ii) filtering, (iii) smoothing.

In a broad sense, both smoothing and filtering are time-dependent DA methods as they both combine the advantages of measurements and model outputs. In a narrower sense, only filtering problems count as time-dependent DA methods, as DA only deals with “real-time” measurements instead of historical ones. In the case that the information in the time series is only propagated forward without any backward loop or window, the time-dependent methods are also termed sequential methods. The prediction problems usually serve as the validation of the smoothing and filtering problems. Traditional batch calibration technologies for hydrological models can be categorized as a smoothing problem (McLaughlin 2002). Filtering technologies usually play an important role in improving the prediction from the calibration; hence, they are widely used in real-time modules of operational models (Divac et al. 2009, McMillan et al. 2013).

2.2.3 Applicability of DA methods to hydrological models

Many approaches have been applied in hydrological DA (Houser et al. 2005, Walker and Houser 2005). The most commonly used methods include the variational method, PFs and the KFs. Although traditionally dominant in numerical weather forecasts, variational methods have not been widely used in hydrological DA. Variational methods such as 3D VAR assume the forecast error statistics are isotropic and largely homogeneous with little variation in time (Houtekamer and Mitchell 1998), though the consideration of the time dimension would overwhelmingly increase the computation burden. It is also complicated to develop the adjoint model for the distributed hydrological models (Clark et al. 2008). Successful cases using variational methods in hydrological DA are mostly based on simpler lumped models (Seo et al. 2003, 2009, Abaza et al. 2014b). Abaza et al. (2014b) compared the classic variational method with EnKF, and they found that the latter is more stable in streamflow assimilation.

PFs are also widely used in hydrology DA (Pham 2001, Moradkhani et al. 2005a, Weerts and El Serafy 2006). One of the major advantages of PFs is that the system does not need to be Gaussian (Liu et al. 2012b). PFs are also better at handling model nonlinearities compared to other sequential methods (Moradkhani et al. 2005a). The fact that PFs use full prior density means this method is more computationally intensive (Weerts and El Serafy 2006). The operational applications of PFs in distributed hydrological models are also limited due to the set-up of the particle numbers (Liu et al. 2012b).

Despite the thriving development of DA methods in “interpolation”, “smoothing” or “filtering”, not many of them are actually extensively validated in “forecast” mode, not to mention in operational hydrological forecasts (Liu et al. 2012b). One of the problems that hydrological DA faces is the short “efficient period”. The updating of model structures and storages generally has a major impact on the forecast only within shorter lead times (El Serafy and Mynett 2004, Knight and Shamseldin 2006). A possible reason for this is that DA deals with uncertainties from different data sources in a statistical way, rather than in a physical or mechanical way (McLaughlin 2002). It somehow “detours” the imperfect model structure and parameters in hydrological models by adding more weight to the model uncertainties, thereby emphasizing observations in the final results.

Numerical Weather Prediction (NWP) is fundamentally an initial problem that is very sensitive to state variation versus the hydrological models, which lean toward a process problem that relies on model structures and parameters. For this reason, DA might not work as well in hydrology as in NWP. However, because of the “conceptual” nature of most hydrological models, DA does have the potential to improve the model forecast either by updating the initial condition or modifying the model parameters.

A potential development direction of hydrological DA is the combination of objective analysis methods and sequential methods. 4D VAR expands the strength of 3D VAR by considering the temporal evolution of variables within the fixed time window. However, depending on the model structure, it could be expensive to calculate the gradients of the cost function (Bin and Ying 2005). The methods that combine the merits of both objective analysis and sequential methods, such as 4D VAR and EnKF, are generally regarded as the most promising DA technology in NWP (Lorenc 2003, Kalnay et al. 2007). However, as far as the authors are aware, such a dominant DA method still does not exist in operational hydrological simulation and forecasting.

3 EKF and EnKF

The linear KF is a classic sequential method. Together with its multiple nonlinear variants, such as the EKF and EnKF, the KF has become a very promising method cluster used in hydrological DA.

3.1 Linear KF

The KF is an optimal estimator that recursively couples the most recent measurements into the linear model to update the model state output (Kalman 1960). Under the assumption that the linear system is a stochastic process with Gaussian noise, the KF produces the best estimation with minimum mean square error.

The KF works on a stochastic system in the form of state space equations (Hamilton 1994):

(1)

(2)

Equation (1) is the model function that propagates state x from step k to step k + 1. x_k is the a posteriori state vector and x_k+1 is the a priori state. M_k is the model function (also known as the dynamic function or dynamic operator). B_k is a linear matrix to convert the dimension of residual vector u_k to state vector x_k. u_k is the “forcing” term of the model in the form of linear residuals. ŋ_k is the model estimation error (covariance matrix ). Sometimes, to simplify the problem, u_k is regarded as part of M_k and its estimation error is included in ŋ_k, and in this case B_k = 0.

Equation (2) is the observation function that relates the state vector x_k to the observation , with observation error ε_k (covariance matrix ). To guarantee the linearity of the system, the observation operator H_k should be linear too. Normally, the number of observation variables is limited while the state variable combinations could be infinite. It is probable that the observation vector dimension (degree of freedom) is much smaller than that of the state vector.

A typical implementation of the KF is described below (Drécourt 2003). At time step 0, set up the initial estimation error covariance matrix of the initial state . The model propagation without an explicit form for u_k is given by:

(3)

(4)

The superscript “a” represents a posteriori estimation and the superscript “f” means a priori estimation.

The state is then updated via:

(5)

where is the innovation and the Kalman gain is calculated with:

(6)

The estimation error is finally updated as:

(7)

With the updated estimation error, the filter can restart from Equations (3) and (4) to begin another loop recursively.

Equation (5) shows that the updated state is a linear combination of the observation and the model estimate (Drécourt 2003). The gain would change between 0 and H⁻¹ depending on the uncertainty comparison between the model outputs and the observations. For example, if the dynamic model is 100% accurate (which is unrealistic), then the new observation would not be considered at all in the updated state, as the gain is 0, and vice versa.

The linear KF tends to be more commonly used in stochastic models (Bolzern et al. 1980, Bergman and Delleur 1985, Szöllősi-Nagy and Mekis 1988) and channel routing problems (Huang 1999), which are easier to linearize (Georgakakos and Bras 1982, Fan 1991, Wu et al. 2008, Sun et al. 2013). However, there are a few successful cases in which the linear KF was applied in complex RR models. Liu et al. (2011) and Lee and Singh (1999) coupled a tank model with a linear KF to estimate the model parameter and outflow respectively. Kim et al. (2005) applied the KF in a one-dimensional physically based distributed model CDRMV3 with the storage amount of the whole watershed as state and outlet discharge as measurement, while a discharge–storage (Q-S) curve was used as the observation equation, in addition to a Monte Carlo ensemble being used as the dynamic equation. The drawback of such operations is that it could be difficult to find the linear relationship between discharge and storage, even at local scale.

3.2 Nonlinear KF

The KF is only legitimate in a linear system where both the model function M_k and the observation function H_k are linear. In the case of a nonlinear system (Drécourt 2003):

(8)

(9)

where is the model function and is the observation function, the linear KF is not applicable even if only one of them is nonlinear.

The hydrological system is a highly complicated nonlinear system with gigantic dimensions. In cases such as stochastic black-box streamflow prediction models, a linear KF might apply. However, it is not feasible to count on a linear KF in DA of distributed/semi-distributed hydrological models.

To cope with the nonlinear problem, many different versions of modified KFs have been developed (Shamir et al. 2010, Muluye 2011, Dumedah and Coulibaly 2012, Chen et al. 2013, de Rosnay et al. 2013, Gharamti and Hoteit 2014, Rakovec et al. 2015). A large portion of these booming new methods or algorithms are based on the EKF and EnKF.

3.2.1 EKF

3.2.1.1 Theory of the EKF

A prerequisite to applying the EKF is that the nonlinear system should be continuously differentiable. Under this condition, Taylor extension is applied at the estimated point. It is therefore possible to obtain the converted linear dynamic matrix by expanding the nonlinear functions at the estimated point. Apply this to both model function and observation function to get Equations (10) and (11) (Walker and Houser 2005):

(10)

(11)

The linearized and are then used in the calculation of error covariance propagation; all subsequent steps are the same as in the linear KF (Lewis et al. 2008). The EKF may be applied to the latest observation back and forth as many times as needed to reduce the linearization errors. However, such recursive operations do not guarantee a better estimation than the one-time EKF, while the computation burden is increased greatly (Puente and Bras 1987).

In many cases, nonlinear systems do not have explicit analytical solutions and the derivatives can only be calculated numerically. As the state is a vector with multiple variables, the linearized matrix can be expressed as a Jacobian matrix of the partial differential functions.

The EKF keeps the first order of the Taylor expansion, while it ignores the higher-order terms. Theoretically, a high-order EKF that keeps more than the first-order terms is more precise, as more information is kept (Tanizaki 1996, Simon 2006). However, it would still be biased because higher-orders terms are still ignored (Tanizaki 1996). In reality, there are rare cases supporting the superiority of higher-order EKF over first-order EKF (Sadeghi and Moshiri 2007, Ermolaev and Volynsky 2014). This is especially true when the analytical solutions of the system are unavailable, because the simultaneous numerical solution of the Jacobian and Hessian matrix may introduce more errors (Roth and Gustafsson 2011, Vittaldev et al. 2012).

3.2.1.2 Application of the EKF

The EKF has been widely used in operational soil moisture analysis by various weather agencies (Hess 2001, de Rosnay et al. 2013, Fairbairn et al. 2014). de Rosnay et al. (2013) described the operational implementation of a simplified EKF (SEKF; Draper et al. 2009, Mahfouf et al. 2009), in which the EKF has been simplified by setting constant background errors in the ECMWF land surface analysis system. The SEKF linearizes the observation operator with finite differences by adding individual perturbations to each element of the model state vector. Such a feature significantly increases the computation cost with the increase of dynamic model resolution. To solve this problem, one should either modify the model structure or calculate the Jacobian matrix from an off-line land surface model simulation (Mahfouf et al. 2009). The comparison between the SEKF and EnKF shows that the SEKF gives similar performance to EnKF; the latter is unable to improve on the former (Fairbairn et al. 2014).

Ge (1984) applied the EKF on a three-layer storage excess runoff model. Depending on the existence of precipitation and the fulfilment of each soil layer, they set 10 different runoff-generation scenarios. To avoid defining a threshold function, each scenario was assigned a state space function. By doing so, it is found that only two state space functions need to be linearized. Taylor expansion was then used with these two functions analytically. The constants from truncation of the Taylor expansion are treated as the linear driven terms of the model function and linear residuals of the observation function, respectively. Neither of them is involved in the state error propagation explicitly, but they are the integral components of the model and observation functions and cannot be neglected.

Kitanidis and Bras (1980a) linearized the Sacramento model using statistical linearization and a “describing function” technique. Instead of linearizing the nonlinear formulas directly, a group of simple linear functions that produce similar outputs with the same inputs as the nonlinear formulas are applied. The focus is thus shifted to the determination of the coefficients of the linear functions using statistical methods. This scheme is somehow “subjective”, because the candidate linear functions are not uniquely dependent on the selection criteria, or the evaluation functions to assess the fitness of the linear functions. Another challenge for this scheme is that it assumes that the inputs follow a constrained Gaussian distribution and, most importantly, they should fall in a very narrow range for the linear function replacement. The computation burden of statistical linearization is also a big issue.

Instead of implementing analytical linearization, Walker and Houser (2001) numerically derived the dynamic transition matrix in the propagation of the estimation error covariance. The state vector contains three intermediate soil moisture parameters and the forecast equations of state were linearized with a first-order Taylor series expansion. The observation chosen is the remote-sensing surface soil moisture, which is related to the state vector through a complicated nonlinear equation that is also linearized with a Taylor expansion. The numerical solution significantly increases the computation cost, but not as much as statistical linearization (Kitanidis and Bras 1980b). Therefore, it is important to provide a state vector of a manageable size. The problem with the numerical solution is the determination of the perturbation of the independent variables: small perturbations may lead to numerical problems, while large perturbations will cause a much greater loss of accuracy, and could also risk hitting the nonlinear threshold (Reichle et al. 2002b).

3.2.1.3 Assessment of EKF

The robustness issue of the EKF is well discussed in hydrological DA, and much of this concerns the divergence problem (Ljung 1979). The divergence of the EKF can be due to a number of reasons (Fitzgerald 1971): The first factor concerns the inappropriate estimation of the model and observation error covariance matrix, which will be discussed extensively in Section 4.2. Secondly, divergence can also be caused by the incorrect estimation of the initial estimation error P. Since the estimation of P at time k relies on P at time k − 1, a large error in initial P could recursively propagate until there is a divergence of the EKF (Tao et al. 2005). The third reason concerns the truncation errors arising from the Taylor expansion: although first-order local linearization is adequate to account for differential nonlinearity in some models (Reichle et al. 2002b), the existence of such errors makes it harder for the model errors to meet the whiteness assumption. Lastly, the errors from the numerical calculation of the Jacobian matrix: the perturbation scale of the independent variables involves subjective trials.

Many efforts have been made to alleviate the divergence of the EKF. Adaptive algorithms (Jwo and Wang 2007) are widely used to prevent the divergence caused by inappropriate error quantifications. Ljung (1979) introduced an innovation model by adding a parameterized Kalman gain term to the original EKF state space equations; he argued that the main reason for the divergence is the lack of coupling between Kalman gain K and model parameters. Tao et al. (2005) rescale the error covariance matrix by multiplying it by a factor; once the divergence is detected, the Kalman gain is frozen until stability is restored.

Compared to other nonlinear methods, the implementation of the EKF is direct and straightforward. It is arguably the de facto standard method to expand the KF to nonlinear systems for its outstanding performance in practice (Rabier and Liu 2003). However, it is not clear if the EKF can outperform other DA methods in nonlinear distributed hydrological models when considering the computational burden of linearization, even though strategies such as low-rank SEEK (singular evolutive extended Kalman; Tuan Pham et al. 1998) can reduce the computational requirements of the EKF.

3.2.2 EnKF

3.2.2.1 Theory of the EnKF

The EnKF is a Monte Carlo approach for nonlinear filtering problems (Evensen 1994a). There are two types of EnKF (Kalnay et al. 2007): perturbed observation EnKF (Evensen 1994b, Burgers et al. 1998) and square-root EnKF (Anderson 2001, Bishop et al. 2001, Whitaker and Hamill 2002). The perturbation of the observations causes additional sampling errors (Kalnay et al. 2007), but this type of EnKF could handle nonlinearities better than the square-root EnKF (Lawson and Hansen 2004).

The EnKF is based on approximation of the conditional probability densities of the state error covariance by a finite large number of randomly generated model trajectories. The EnKF does not need any derivation of the model operator or observation operator as performed in the EKF. Instead, it generates a set of realizations (the ensemble) and propagates them through the model operator independently. Then, it derives the a priori state error covariance through the statistical analysis of the ensemble (Evensen 1994b).

The general idea of a perturbed observation EnKF is demonstrated in Figure 2. In the forecast step, the a priori state ensemble is created by adding random perturbations to the best estimate of the initial state. Then the a priori error covariance matrices are approximated from the state and model output ensemble error matrices (Gillijns et al. 2006). In the analysis step, the state ensemble x are integrated forward in parallel, based on the original model. The observation is perturbed by adding Gaussian noise of measurement errors (Burgers et al. 1998). Then, the a posteriori error covariance matrix can be estimated from the ensemble error matrix of the a posteriori states. Meanwhile, the Kalman gain is calculated from the forecast error covariance matrices (Reichle et al. 2002b, Moradkhani et al. 2005b).

Figure 2. Schematic of the EnKF with perturbed observations (Moradkhani et al. 2005b). Reprinted from Dual state–parameter estimation of hydrological models using ensemble Kalman filter. H. Moradkhani, S. Sorooshian, H. V. Gupta, P. R. Houser. Advances in Water Resources 28 (2005), 135–147. Copyright (2005), with permission from Elsevier.

3.2.2.2 Application of the EnKF

Although the propagation of state error covariance does not require the linearization of the model, the parallel computation of the ensemble means the EnKF is an expensive method to use. Nevertheless, the EnKF is a widespread approach in hydrological DA.

Weerts and El Serafy (2006) compared the EnKF with two PFs in a flood forecast with the conceptual RR model HBV-96, and found that the EnKF outperforms both PFs. Dumedah and Coulibaly (2013) compared evolutionary data assimilation (EDA) with methods based on the integration of Pareto optimality into both an EnKF (ParetoEnKF) and a PF (ParetoPF). Rafieeinasab et al. (2014) compared the EnKF with a maximum likelihood ensemble filter, which is an ensemble extension of variational assimilation. It is argued that the former is less sensitive to observation, model errors, and parameter uncertainties, while it does perform reasonably well with a smaller ensemble size. Borup et al. (2015) discussed a procedure to utilize the EnKF in the assimilation of observations that are “out-of-range”, and defined this method as “partial EnKF”. The positive results indicate that the EnKF is versatile enough to take advantage of the imperfect but precious observation information, instead of simply wasting it.

Madsen et al. (2003) presented a procedure to assimilate observed water levels and fluxes in the MIKE 11 Flood Forecasting system using the EnKF. Borup (2014) demonstrated the application of the EnKF in assimilating water-level and flow observations into distributed urban drainage models. Liu et al. (2015) used the EnKF to assimilate AMSR-E snow depth into a land surface model for streamflow predictions. Xie and Zhang (2010) implemented synthetic simulation experiments for application of the EnKF with a Soil and Water Assessment Tool (SWAT) model. Sensitivity analysis with regard to error specification, initial realization and ensemble size were also demonstrated in their paper. Applications of the EnKF in Soil and Water Assessment Tool (SWAT) were also extensively studied in other works (Han et al. 2012, Lei et al. 2014).

Clark et al. (2008) attempted to apply the EnKF to assimilate streamflow observations into a distributed hydrological model; however, the researchers found it unsuccessful due to the inappropriate selection of state variables. Reichle et al. (2002a) assimilated remote-sensing soil moisture into a land surface model with a very large state vector and a small ensemble number. They argued that a large ensemble number is needed to obtain a robust error variance estimation. Shi et al. (2014) found that the EnKF assimilation of multivariate observations applies strong constraints to parameter estimation in a physically based land surface hydrological model. Abaza et al. (2014a) introduced hyper-parameter perturbation factors by comparing the H-EPS spread to its mean forecast error for perturbation of the system in the EnKF. This research outlined the importance of input replicate generation in the implementation of the EnKF. Panzeri et al. (2015) described the implementation of stochastic moment equations (MEs) based on the EnKF in groundwater flow simulation; the researchers found it more efficient than the traditional EnKF.

3.2.2.3 Assessment of the EnKF

Compared to the EKF, the EnKF is less likely to diverge because the Taylor expansion and Jacobian matrix calculation are averted. However, the EnKF accounts for more model errors, as there is no truncation of higher-order terms as in Taylor expansion (Reichle et al. 2002b). Considering the computational burden of the linearization, it is unrealistic to put all the representative variables from the calculation unit of the distributed model into the state vector (Gillijns et al. 2006). Therefore, it is difficult to consider both the horizontal correlations in the model, and the measurement errors in the EKF. However, this is not an issue for the EnKF (Reichle et al. 2002b).

Computational efficiency might be the primary concern before one applies the EnKF. A major drawback of the EnKF is that it assumes that the prior PDFs of the model states are Gaussian; hence, the posterior states are only determined by the first two moments of the prior density (Weerts and El Serafy 2006). Theoretically, the larger the ensemble number, the more accurate the estimation of the mean and variance of the prior density (Madsen et al. 2003, Rasmussen et al. 2015). It is argued that the EnKF estimate would not converge to an optimum unless the ensemble size is large, especially in the case of high-dimensional problems (Reichle et al. 2002a, 2002b, Crow and Wood 2003). One of the preconditions to consider the EnKF superior to the EKF is that the ensemble size should reach a certain level, so that the major errors are statistical noise, rather than a closure problem or unbounded error variance growth encountered by the EKF (Evensen 1997). In reality, parallel running of distributed hydrological models is a huge challenge to computational capacity. Although the determination of reasonable ensemble size is regarded as a case by case problem (Xie and Zhang 2010), more efforts are required to derive efficient algorithms to trade between ensemble size and balanced forecast analyses, as well as representative error statistics (Mitchell et al. 2002).

4 Implementation of KFs

4.1 Filter elements selection

An incomplete list of recent hydrological DA works has been made to summarize their features, including the chosen DA methods, hydrological models, state, parameters (if any), and observations (see Table 1). Depending on the assimilation purpose and the available observation, the selection of which state variables to assimilate can be very flexible. A state variable can be either intermediate (e.g. soil moisture content, snow–water equivalent) or prognostic (e.g. runoff/discharge). Soil moisture is by far the most popular state variable, because of its crucial role in describing the model with state space equations. Nowadays, there are many remote-sensing products that can provide temporally and spatially continuous soil moisture observations (Moradkhani 2008). On the other hand, streamflow forecasting is one of the ultimate goals of hydrological modelling. With the easy availability of the streamflow observations, it is also a widely used state variable.

Table 1. Literature summary of state and parameter(s) selection.

Reference	Method	Model	State	Parameter(s)	Observations
Reichle et al. (2002a)	EnKF	Reichle et al. (2001)	Soil moisture and soil temperature	–	L-band microwave radio brightness
Madsen et al. (2003)	EnKF	MIKE 11	Water level and flux	–	Water level and flux
Moradkhani et al. (2005b)	EnKF	HyMOD	Water storage	Rq, Rs, α, bexp, Cmax	Streamflow
Vrugt et al. (2006)	EnKF	SAC-SMA	Six soil-zone reservoir variables, three routing variables	14 time-invariant parameters	Streamflow
Weerts and El Serafy (2006)	EnKF, PF	HBV96	Soil moisture, snow water, interception storage, lower and upper zone, etc.	–	Discharge
Clark et al. (2008)	EnKF	Topnet	Soil storage, aquifer storage, surface storage, streamflow prediction	–	Streamflow
Xie and Zhang (2010)	EnKF	SWAT	Runoff, water storage, evapotranspiration, etc.	CN2	Streamflow, soil water content, evapotranspiration
Dechant and Moradkhani (2011)	EnKF, PF	Snow 17, SAC-SMA	SWE	All model parameters	Microwave radiance
Chen et al. (2011)	EnKF	SWAT	Soil moisture	–	Surface soil moisture
Cammalleri and Ciraolo (2012)	EnKF	SVAT	Soil moisture	All model parameters	Surface temperature, latent heat flux, soil water content
Han et al. (2012)	EnKF	SWAT	Soil moisture	–	Surface soil moisture
Trudel et al. (2014)	EnKF	CATHY	Pressure head, incoming and outgoing discharge	–	Streamflow
Samuel et al. (2014)	EnKF	SAC-SMA	Soil moisture content, SWE, etc.	10 time-variant parameters	Streamflow, soil moisture
Kitanidis and Bras (1980b)	EKF	NWS	Soil water content, additional impervious storage, channel storage, routing coefficient	–	Streamflow
Puente and Bras (1987)	EKF	Conceptual lumped model	Water equivalent mass in cloud, soil water content, channel reservoir storage, etc.	–	Discharge, precipitation
Francois et al. (2003)	EKF	GRKAL	Soil moisture	–	Surface soil moisture
Aubert et al. (2003)	EKF	GR4J	Water level in soil reservoir and water level in routing reservoir	–	Soil moisture, streamflow
Kumar and Kaleita (2003)	EKF	LSM	Soil temperature	–	Near-surface soil temperature
Jonsdottir et al. (2006)	EKF, maximum likelihood	Stochastic model similar to HBV and NAM	Temperature, snow, upper and lower water content	All model parameters	Discharge
Lü et al. (2011)	EKF, PSO	1-D Richards’ equation	Soil moisture content	Ks, θr, θs, α, n	Surface soil moisture

PSO, particle swarm optimization.

In most cases, the model parameters are assumed to be time invariant once they are calibrated. However, it is possible to treat some parameters as time variant, and update them together with the state variables (Samuel et al. 2014), even though this is still not as widespread as it should be. The selection of the parameters is less flexible than that of the state variables because it should coordinate with the state variable update. Also, the parameters selected should be sensitive enough to reflect the update (Xie and Zhang 2013). In the case of a model being simple enough, it is possible to include all the parameters in an assimilation; however, this would present a huge challenge for complex models. He et al. (2012) developed an Integrated Sensitivity and Uncertainty analysis Framework (ISURF) to screen and identify the sensitive model parameters and assess the uncertainty structure of model parameters in the EnKF.

4.2 Determination of errors

Another key issue in the implementation of the KF is the determination of model and observation errors. The model errors are more difficult to describe as compared to the observation errors, because the latter can sometimes be predefined based on the measurement, sampling methods or empirical formulas (McMillan et al. 2010). The overestimation of model errors reduces confidence in the model, and thus the filter would overly rely on observations; to the contrary, the underestimation of model errors exaggerates the accuracy of the model and also wastes information from new observations (Kitanidis and Bras 1980b). The improper selection of model errors may lead to unacceptable results or a divergence of the KF. Puente and Bras (1987) argued that proper error quantification of the model is even more important than the selection of the DA methods.

Clark et al. (2008) quantified the errors in streamflow measurements as a fixed proportion of the discharge observations, while the model errors were generated by perturbing the precipitation forcing and model states (e.g. soil moisture and aquifer storage) with temporal varying perturbations. However, the more general method to realize the time variant estimation of model and observation error is by employing adaptive filtering. If properly used, the adaptive filtering may partly mitigate the influence of the inaccurate (if not incorrect) set-up of the initial model and observation errors. The adaptive filters are divided into four categories: Bayesian, maximum likely hood, correlation and covariance matching (Mehra 1972).

Bayesian methods aim to obtain recursive equations for the a posteriori probability density of states and errors given the new observations. The Bayesian approach is considered as the theoretical basis of all other approaches (Sarkka and Nummenmaa 2009). Examples of Bayesian approaches include state augmentation methods (Wan and Van Der Merwe 2000) and interacting multiple models methods (Rong Li and Bar-Shalom 1994). To improve the operational flood forecast, Li et al. (2014) used a maximum a posteriori method constrained by prior information drawn from flow-gauging data to estimate model and observation errors in the Ensemble Kalman Smoother (EnKS) assimilation of discharge.

The Sage-Husa method (Sage and Husa 1969, Sage and Wakefield 1972), which is generally used in linear discrete systems, is a maximum posterior probability method that recursively estimates the mean and covariance of the unknown system and observation noise. The Sage-Husa method is based on the assumption that the system is stationary; therefore, the model and observation error covariance will eventually converge to constants. In the case of non-stationary systems, Duan et al. (2003) proposed to rescale the model errors with a newly defined time-varying factor to avoid the likely divergence of the filter.

The EKF and EnKF are both suboptimal filters, this is mainly due to the representativeness of errors (Mitchell et al. 2002, Reichle et al. 2002b). For suboptimal filters with non-white innovations, Mehra (1970) suggested an auto-correlation function of the innovation process to obtain asymptotically unbiased and consistent estimates of model and observation errors. However, this method requires the number of the unknown elements in the model error to be no greater than the state vector size by observation vector size; in addition, the structure of the model error should also be predefined. This method does not guarantee the uniqueness of the estimated model error under the given conditions (Odelson et al. 2006). Thus, the representativeness of the sample correlation function is questioned, and the minimum covariance may not be reached due to the likely correlation between the observations (Neethling and Young 1974).

The basic idea behind the covariance matching method is to make the residuals consistent with their theoretical covariances (Mehra 1972). This method is particularly popular in atmospheric data assimilations (Dee 1995). Reichle et al. (2008) modified the covariance matching approach and used it in surface soil moisture assimilation. In a nutshell, their modified method assumes that the internal diagnostics of the assimilation system are consistent with the values expected from the input parameters; as a result, the adaptiveness could be realized by adjusting the “adaptive tuning factor”, which is the ratio of the true input error variance and their initial values. Like most other covariance matching methods (Mehra 1970), this method works better for the identification of observation errors than model errors. It is also believed that covariance matching techniques give biased estimates of the true covariance (Odelson et al. 2006).

In general, adaptive update methods normally work well in linear systems with Gaussian errors, under stationary conditions (Reichle et al. 2008). However, some of them exploit the “optimality” of the KF, especially in the assumption of the whiteness in the innovation sequences (Kitanidis and Bras 1980b). In the case of non-stationary systems or non-white errors, computation becomes extremely tedious and complicated. Crow and Reichle (2008) compared four adaptive filtering schemes in land data assimilation and subsequently found they all have the problem of low convergence. It is also argued that the current self-adjustment or adaptive methods are feasible only when system dimensionality is reduced (McLaughlin 2002). To seek a parsimonious DA approach that requires a minimum assumption to describe the unknown model and observation error characteristics, Vrugt et al. (2005) and Crow and Yilmaz (2014) developed the Auto-Tuned Land Data Assimilation System (ATLAS) for the integration of two remote-sensing soil moisture products into a water balance model. They combined four separate adaptive filtering solutions (the triple collocation, innovation, merged, and red modelling error solutions) to estimate the observation and model errors. The application of ATLAS to a simple forecast model leads to an improved surface soil moisture analysis, although its power with more complex models is yet to be investigated (Crow and Yilmaz 2014).

4.3 Parameter estimation

Due to limited knowledge of hydrological processes, it is impossible to predetermine all the parameters in hydrological models. The most common practice first involves the initialization of the parameters, which are then adjusted with batch calibration until certain criteria are met. Once determined, the calibrated parameters are assumed to be consistent in future simulations and are rarely updated again.

The KF is designed to update the state of the system for a better initialization of the forecast. However, it is not easy to find a state vector with a significant representativeness to update the overall estimation. Moreover, a pure state update neglects the propagation of the uncertainties in model parameters (Lü et al. 2011). Thus, in order to simultaneously update the model parameter with the state, it is desirable to account for both the uncertainties in parameters, and also state estimation (Liu and Gupta 2007). There are three main schemes to estimate the model parameters: state augmentation, dual state-parameter estimation and hybrid solutions.

4.3.1 State augmentation

State augmentation technology is realized by adding parameters to the original state vector; thus, both the original dynamic model and observation model would be modified (Hendricks Franssen and Kinzelbach 2008). It is rare that the relationship between state variables and model parameters are linear; therefore, the linear KF is normally inapplicable when the state vector is augmented. In the case of augmented EnKF, Hendricks Franssen and Kinzelbach (2008) argued that the state and parameters can be jointly updated with either an iterative or a non-iterative approach. In the case of the EKF, one can normally implement the Taylor expansion and Jacobian matrix methods as usual. However, it is worth pointing out that the parameters are usually independent from prognostic state variables, which means the corresponding derivatives in the Jacobian matrix are zeros. On the other hand, the prognostic state variables are not independent of the parameters.

Based on Equations (8) and (9), supposing is the augmented new “state vector” that contains both the original state vector and parameters , the augmentation process can be demonstrated as follows (Wang and Wang 1985, Gharamti and Hoteit 2014):

(12)

where is the new model operator and is the new model error. Assuming the parameters are constants, the new model error can be expressed as:

(13)

where

(14)

where is the original model error covariance and 0 is zero matrix.

Meanwhile, the observation operator is also transformed to:

(15)

where is the original observation function.

Equations (13) and (14) are based on the assumption that the parameters are constants. In the case of varying parameters, the following equations should be used (Wang et al. 2009):

(16)

where is the parameter error. Assuming is independent of , then:

(17)

where is the covariance matrix of the parameter errors.

The augmentation of parameters could lead to more complex numerical solutions to the Jacobian matrix. When numerical techniques such as forward, backward or central Euler schemes are utilized (Chang and Latif 2009), the EKF might be very sensitive to the chosen step size of the parameters. The EKF is also sensitive to truncation errors and round-off errors of the parameters, which therefore makes it more difficult to choose the optimal step size (Gharamti and Hoteit 2014). Although the EnKF is more capable of handling state vectors with a large degree of freedom, the increase in the degree of freedom caused by parameter augmentation, together with the spurious long-range correlations between state variables and parameters, may contribute to false parameter estimation and offset the benefits of parameter updates in the EnKF assimilation (Reichle and Koster 2003, Xie and Zhang 2013). Reichle and Koster (2003) applied localization (Keppenne and Rienecker 2002) to constrain the correlations of the state vector elements beyond a certain separation distance. Xie and Zhang (2013) introduced a partitioned parameter update scheme, which is logically similar to the implementation of partially derivative repeated updates on each parameter segment, while the rest remain fixed. This new scheme is believed to be superior to regular state augmentation (Xie and Zhang 2013, Rakovec et al. 2015). Two major drawbacks arise for this method: first, the repetitive run of the model not only significantly increases the computation demand, but may also accumulate model errors; second, it is found the update order of the parameter segments may influence the assimilation results (Xie and Zhang 2013).

4.3.2 Dual state-parameter filter

The dual state-parameter filter method dates back to the development of the Mutually Interactive State-Parameter (MISP) method (Todini 1978). In this method, two KFs are run simultaneously, one for the state update and one for the model parameter update. For each time step, the two filters are run alternately once the new observation is available: the update of the state benefits from the update of the parameters through the update of the model function, and the update of the parameters benefits from the update of the state from the update of the innovation (Bergman and Delleur 1985, Wang and Wang 1985).

The equations of a dual state-parameter EKF with a nonlinear model function and a nonlinear observation function are provided, based on the following calculation (Haykin 2001):

First, the parameters and parameter estimate errors are propagated:

(18)

(19)

where is the a priori parameters estimate, is the a posteriori parameters estimate, is the a priori estimate error of the parameters, is the a posteriori estimate error of the parameters, and is a “forgetting factor” (Haykin 2001, Ji and Brown 2009).

Then, after running the state filter:

(20)

(21)

(22)

(23)

(24)

where is the a priori estimate error of the state, is the a posteriori estimate error of the state, is the Kalman gain for the state. and are the linearized model and observation functions calculated with Equation (10) and (11). The rest of the variables have the same meanings as shown in Equations (3)–(7). A major difference between the dual state-parameter filters and the state only filter is that the model function in the former is updated with the parameters.

Define:

(25)

(26)

The parameters can be updated:

(27)

(28)

The application of the dual state-parameter method is less common compared to state augmentation. Moradkhani et al. (2005b) described the application of dual state-parameter filters in the EnKF. In that case, the Kalman gain of the state is calculated with the cross-covariance of the states ensemble and the prediction ensemble, instead of using Equation (22). The Kalman gain of the parameters is calculated through the cross-covariance of the parameter ensemble and prediction ensemble instead of using Equation (27) (Moradkhani et al. 2005b).

4.3.3 Hybrid solutions

The last category of the parameter update is named “hybrid solutions”, because it involves a third party scheme other than KFs. Lü et al. (2011) coupled the EKF with the optimal parameter estimation, which is realized with a PSO, to obtain the dual state-parameter estimation of root-zone soil moisture within Richards’ equation. This coupling was found to indeed improve the simulation results. Vrugt et al. (2005, 2006) presented the Simultaneous Optimization and Data Assimilation method (SODA), which uses the EnKF to recursively update model states conditioned on an assumed parameter set within the inner loop, while additionally estimating time-invariant values in an outer global optimization loop using the shuffled complex evolution metropolis stochastic-ensemble optimization approach. It is proved that SODA not only improves the estimate of model parameters and state variables, but also creates reliable model prediction uncertainty bounds and a time series of valuable state and output innovations.

4.4 Routing problem

The EKF and EnKF are designed to propagate forward, with the model’s forecasted state predictions being updated at the same time as the observation is obtained (Crow and Ryu 2009). In reality, the gap between the continuous hydrological states and the discrete observation usually makes this assumption hard to satisfy. For example, when the hourly RR model is used in DA, the discrete discharge “observation” at the watershed outlet is controlled by not only the watershed water storage “state” at the current hour, but also that of the last few hours, when considering the lag time that the water within the watershed took to travel to the outlet (Pauwels and De Lannoy 2009). Neglecting the lag time may cause serious consequences if a distributed model is used in a large watershed, while the underground water flow is also taken into account. This is an issue that all sequential methods face when they are used in hydrological DA. However, it is mostly discussed in the context of the EnKF for its widespread application in such problems.

Weerts and El Serafy (2006) skipped the routing problem by shifting the state forward for a period, during which time most of the rainfall leaves the watershed to reach the outlet. In other words, the future observation is “borrowed” to be used in the current assimilation. As simple as it is, this method neglects the baseflow component in the channel; also, the determination of the “shift time” could be subjective, depending on the size of the watershed, as it is obviously different from the time of concentration. Pauwels and De Lannoy (2006, 2009) suggested a retrospective EnKF in which the current state vector is augmented with a few past states. The number of augmented historical states is determined by the time of concentration. Assuming the assimilation of each time step provides the initial condition for the next step, the EnKF is reapplied forward from the start of the concentration through the current time step. Unfortunately, their results do not seem solid enough to support the efficiency of such methods, due to the accumulation of model errors after each retrospective run (Pauwels and De Lannoy 2009). Meanwhile, McMillan et al. (2013) suggested a recursive EnKF that iteratively updates past and present model states with lags up to the concentration time of the catchment to improve the initial condition, hence the streamflow forecast. This method is found more stable compared to the EnKF, while the computation cost could be an issue if the model is complicated and the watershed is reasonably large. The EnKS uses the EnKF analysis as a first estimate, and propagates the updating backward in time. The state covariance matrix is extended to the past time steps within an analysis time (Evensen and Van Leeuwen 2000, Dunne and Entekhabi 2005). Li et al. (2013, 2014) found that the EnKS is superior to the EnKF only when the soil moisture was updated without routing storage being considered. The EnKS is found to be considerably superior to the EnKF when verified in the hypothetical forecasting mode. Sakov et al. (2010) introduced Asynchronous EnKF (AEnKF), which is essentially equivalent to 4D VAR without a tangent linear or adjoint model, but does however require the forward integration of the model (Rakovec et al. 2015). The AEnKF is believed to be more suitable for large-scale systems compared to the EnKS because it does not require multiple ensemble updates. Rakovec et al. (2015) demonstrated the application of the AEnKF by augmenting the state vector with past-forecasted observations, and argued that the AEnKF can be considered an effective method for model state updating from operational aspects.

It is worth pointing out that the prerequisite to use the methods mentioned above is that the routing scheme used in the model is the unit hydrograph, instead of the storage routing scheme (i.e. linear reservoirs). The lag time is not an issue when the storage routing scheme is applied (Moradkhani et al. 2005a, Vrugt et al. 2006).

In conclusion, the mainstream of the existing schemes outlines (1) augmenting the state vector with historical observations and (2) a retrospective run of the EnKF. The second point may slightly violate the conception of “filtering”, but nevertheless it seems it has been well accepted as the rule of thumb to cope with the routing issue.

5 Summary and conclusions

This paper examined and discussed the latest developments of Kalman-type hydrological DA. Major attention has been paid to the implementation and assessment of KF-type DA, especially that of the EKF and EnKF. Nevertheless, many of these issues also apply to other hydrological DA methods, while the intensive research on KFs provides a better platform to discuss them.

The EKF features direct linearization of the system by applying Taylor expansion and solving the Jacobian matrix. The robustness issue of the EKF is amply discussed in hydrological DA, and much of this concerns the divergence problem. Many reasons can cause divergence of the EKF: erroneous determination of the model and observation error, improper initialization, truncation error of Taylor expansion, and numerical calculation errors of the Jacobian matrix, among others. The EnKF is a Monte Carlo method based on a Gaussian distribution assumption of the ensemble. Despite many scholars regarding the EnKF as a solution to the computational issues the EKF faces, the computational efficiency of the EnKF itself is a primary concern. Determination of a reasonable ensemble size is regarded as a case by case problem, although more work is needed on the trade-off between ensemble size and forecast accuracy, as well as the error representativeness.

It is critical to determine the model and observation error in KFs. Adaptive filtering is one of the possible solutions to this issue. Four categories of adaptive filters are summarized: Bayesian, maximum likelihood, correlation and covariance matching. Despite the adaptive filters having already been applied to stationary linear systems with Gaussian errors, their applicability is still questioned with regard to non-stationary systems or systems with non-white errors. It is also believed that adaptive methods only work well when the model dimensions are significantly reduced.

Joint assimilation of state and parameter has been adopted in many hydrological DA studies. There are three major schemes to realize this: state augmentation, dual filter and hybrid methods. In the foreseeable future, it is unlikely that the joint estimation of state and parameter would replace batch calibration completely. Nevertheless, such an operation may help to extend the “efficient period”, considering that the short period of influence is one of the obstacles that limits the application of KF DA in operational hydrological forecasting (Muluye 2011, Agboma and Lye 2014, Grassi and De Magistris 2014).

Despite fruitful achievements in Kalman-type hydrological DA, some challenges remain. These include, but are not limited to:

1. The quantification of model and observation errors. Even though the adaptive filters may alleviate the influence of improper initial error quantifications, it is not clear as to what extent they may alleviate. Most of the current research still tends to quantify model and observation errors using trial and error or empirical equations without further updates to adaptive filters.

2. The assimilation of model parameters through a KF is still in its infancy. The “update” of the parameters somehow contradicts the hypothesis that the parameters are “calibrated” constants, yet there is a lack of specific comparison studies between the batch calibration and the KF assimilation update for the determination of the parameters.

3. The application of Kalman hydrological DA in large-scale watersheds. The major difference between large and small watersheds is the temporal mismatch of the outlet streamflow observation and the update of watershed and channel water storage. The update of watershed and channel water storage may not be verified with the outlet streamflow observation immediately, considering the significant travel time for the water to reach the outlet in a large watershed.

4. The unknown true state. Most of the EnKF applications utilize synthetic simulation results to represent the true state of the system (some even use them as quasi-observations; Crow and Wood 2003). In the forecast, both EKF and EnKF applications eventually compare their assimilation results with the observations, of which the historical parts have already been assimilated. This indicates that the observations are consciously or unconsciously presumed identical to the true state, and are free of errors (Ahsan and O’Connor 1994). The identification of the “true state” is also a key issue to promote the Kalman hydrological DA from the design and experiment stage to the operational stage.

5. The handling of nonlinearity. The EKF and EnKF are two of the most recognized nonlinear variants of the linear KF, yet they are both suboptimal and far from perfection. More methods are desired that could balance the assimilation accuracy and computation cost, as well as the implementation simplicity.

6. The coupling of DA with data-driven models (Solomatine and Ostfeld 2008). The traditional process-based RR models place greater emphasis on model structure and parameters; these hinder the coupling of Kalman hydrological DA in at least in two ways: failure to create an interface to merge multiple source observations (e.g. the remote-sensing products), and the models being less sensitive to the initial conditions compared to input and parameters.

Disclosure statement

No potential conflict of interest was reported by the authors.