The streamflow estimation using the Xinanjiang rainfall runoff model and dual state-parameter estimation method

Journal of Hydrology 480 (2013) 102–114

Contents lists available at SciVerse ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

The streamflow estimation using the Xinanjiang rainfall runoff model and dual state-parameter estimation method Haishen Lü a,b,⇑, Ting Hou a, Robert Horton c, Yonghua Zhu a, Xi Chen a, Yangwen Jia d, Wen Wang a, Xiaolei Fu a a

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, College of Hydrology and Water Resources, Hohai University, Nanjing 210098, China Department of Earth and Environmental Sciences, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada Department of Agronomy, Iowa State University, Ames, IA 50011, USA d Department of Water Resources, Institute of Water Resources and Hydropower Research, 20 Che-Gong-Zhuang West Road, Beijing 100044, China b c

a r t i c l e

i n f o

Article history: Received 8 September 2012 Received in revised form 26 November 2012 Accepted 8 December 2012 Available online 19 December 2012 This manuscript was handled by Konstantine P. Georgakakos, Editor-in-Chief, with the assistance of Eddy Y. Zeng, Associate Editor Keywords: Xinanjiang rainfall runoff model Streamflow estimation Particle swarm optimization Parameter sensitivity Ensemble Kalman filter

s u m m a r y To accurately estimate floods with hydrological models, the model parameters and the initial state variables must be known. Good estimations of parameters and initial state variables are required to enable the models to make accurate estimations. The Xinanjiang rainfall-runoff (RR) model has been widely used in humid and semi-humid regions in China. In this paper, we evaluate the sensitivity of the Xinanjiang RR model parameters and the correlation between the state variables and the output streamflow. In order to reduce the impact of streamflow data error, model structural error and parameter uncertainty, the Xinanjiang RR model is coupled with intelligent optimization algorithms and a data assimilation method to estimate the streamflow in the Luo River in China that was first considered as gauged, and then as ungauged for parameter calibration. Model parameters are estimated in batch using a particle swarm optimization algorithm and three variations of ensemble Kalman Filter data assimilation (the state variable assimilation, parameter assimilation and the dual assimilation at the same time) for the Luo River, and the ungauged basin. In this case, the model parameters are set equal to median values. The results show that when parameter values are determined by an optimization algorithm using 10 years of data, the dual ensemble Kalman Filter notably improves the simulated results. When the parameters adopt the median, the state variable assimilation has little effect on the estimation results, but the parameter assimilation positively affects the simulated results. The dual ensemble Kalman Filter also notably improves the simulated results. The time scale of assimilation has little effect on simulation results for the dual assimilation, but it has a large effect on the state variable assimilation and the parameter assimilation. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction In recent decades, a large number of rainfall-runoff (RR) models has been developed, representing the main physical processes involved in the transformation of rainfall into runoff (Jones et al., 2008; Lindstrom et al., 1997; Moore, 2007; Yu, 2000). Among these models, some of them are lumped conceptual RR model such as the Sacramento model (Burnash et al., 1973), the HBV-96 model (Lindstrom et al., 1997) and the Xinanjiang model (Zhao, 1992). Such models were used for various applications, including hydrology process studies (Bingeman et al., 2006; Li et al., 2009), runoff estimates (Xu, 1999), soil moisture variations (Liu et al., 2007; Lü

⇑ Corresponding author at: Department of Earth and Environmental Sciences, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada. Tel.: +1 519 8884567x37343 (O); fax: +1 519 7467484. E-mail address: [email protected] (H. Lü). 0022-1694/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2012.12.011

et al., 2009; Yu et al., 2010) and statistical analyses of hydrological processes in ungauged basins (Akhtar et al., 2008; Li et al., 2007). Lumped conceptual RR models often include state variable relationships as heuristic or empirical relationships (Kuczera, 1998). These RR models typically have 10 or more parameters, some model input, some interim state variables and some model output (Vrugt et al., 2006). The parameters cannot be directly obtained from the basin characteristics using existing measured streamflow data. Rather, by gradually adjusting the model parameters to match data, a set of accurate model parameters are determined (McCabe et al., 2005; Serbert, 2000). The Xinanjiang RR model has been widely applied in humid and semihumid regions in China (Cheng et al., 2002; Jayawardena and Zhou, 2000; Ju et al., 2009; Li et al., 2009; Zhao, 1992; Zhao et al., 1980). Gan et al. (1997) reported consistent comparisons of the Xinanjiang model with the Pitman model of South Africa, the Sacrament model of USA, the NAM model of Denmark and the SMAR

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

model of Ireland. This model has 14 parameters whose values must be specified. Some of these parameter values can be estimated directly from knowledge of physical watershed characteristics (Ju et al., 2009; Li et al., 2009; Zhao, 1992). However, most parameters represent ‘‘effective’’ watershed properties that cannot, in practice, be measured via direct observations in the field. It is a common practice to estimate parameter values with a historically observed streamflow hydrograph using either an ‘‘expert’’ trial-and-error approach or an automated search algorithm (Thiemann et al., 2001). The parameters estimated in this manner are spatially and temporally lumped representations of heterogeneous watershed properties. Using historical data lacks the flexibility to investigate possible temporal evolution of the model parameter, Thiemann et al. (2001) emphasized another limitation of batch calibration in hydrological prediction of an ungauged watershed when the lack of sufficient historical data makes the batch method infeasible. The results of these factors have led to a number of publications on the treatment of uncertainty caused by model parameters (Clark et al., 2008; Duan et al., 1992, 1993; Misirli et al., 2003; Moradkhani et al., 2005a; Salamon and Feyen, 2009; Thiemann et al., 2001; Xie and Zhang, 2010). Hydrological model studies have focused on the following two points: (1) improving estimation of model parameters using optimization methods; (2) improving the estimation with time of state variables when the model parameters are pre-determined. Common batch calibration techniques only consider parameter uncertainty, while ignoring the uncertainty of the input, output and model structure. All errors (including input variables, output variables, model structure and parameters) are owed to the model parameters. This approach results in a large error between simulated values and observed values. Sequential data assimilation techniques provide a means of explicitly taking account of input, model parameters and output uncertainties (Evensen, 2009). The Kalman filter (Kalman, 1960) (KF), developed for linear systems with Gaussian uncertainties, is an early data assimilation technique. For nonlinear models, KF was later adapted into the extended Kalman filter (EKF). KF and EKF have been widely used in hydrological models (Lü et al., 2010; Moradkhani et al., 2005b; Weerts and Serafy, 2006; Xie and Zhang, 2010). If nonlinearity is strong, EKF is no longer applicable. This led to the development of the ensemble Kalman filter (EnKF) where errors are allowed to evolve with the nonlinear model equations by performing an ensemble of model runs. The EnKF is widely used in RR models (Evensen, 1994; Liu and Gupta, 2007; Vrugt et al., 2006; Weerts and Serafy, 2006; Xie and Zhang, 2010). Weerts and Serafy (2006) indicate that the EnKF is quite robust. The Xinanjiang RR model was developed by Zhao et al. (1980) and Zhao (1992), and was successfully applied in the Southern area of China (Ju et al., 2009). It is important to discover how to determine Xinanjiang model parameters with advanced optimization algorithms. The objectives of this paper are to: (1) investigate the sensitivity of the Xinanjiang model parameters and the correlation between the state variables and the output streamflow; (2) determine the model parameters using a particle swarm optimization (PSO) algorithm; (3) couple the Xinanjiang RR model and the EnKF data assimilation method; (4) use the coupled model at an ungauged basin.

2. Methodology 2.1. Study area description Site 1: The main study area (see Fig. 1) is the Luo River basin located in Guangdong, China. Some hydrology and climate

103

statistical characteristics are summarized in Table 1. Four rain gauges and evaporation stations were selected in the main streams of the basins: Nangao, Luoyan, Dingyang and Nanwan. Precipitation and evaporation data from the four rain gauges and from four evaporation stations in the area were used for daily runoff simulations. Daily rainfall and evaporation data are available from these four stations and daily streamflow data are available from the Nangao hydrologic control stations for the 16-year period of 1988– 2003. Type E-601 evaporation pans were used to observe the amount of evaporation during this period. Site 2: The Shiguan River basin (Fig. 1) is located in the southern part of the Huaihe River Basin which is a closed basin. Some hydrology and climate statistical characteristics are summarized in Table 1. Jiangji is the basin outlet. The basin slopes from south to north with the elevation changing from 1536 to 18.7 m over a horizontal extent of about 128 km. The northern part is characterized by hills, hillocks and river valley plain. The southern part is mainly plain terrain. The Shiguan River is an important tributary of the Huaihe River and originates from the northern slope of the Dabie Mountain in the south Huaihe River Basin. The river system has two major tributaries: the Shi River to the east and the Guan River to the west, and they merge near the Jiangji outlet (Fig. 1). The two tributaries are regulated by the Meishan and Nianyushan reservoirs over watershed areas of 1970 and 924 km2, respectively, with corresponding travel distances of 88 and 89 km from reservoir outflow to the Jiangji outlet. Climatologically, it lies in the warm semi-humid monsoon region. The sub-basin exhibits a diversified natural topography, land surface cover, soil texture, hydrogeology, hydrology and climatology. Precipitation mainly occurs in the period from mid-May to mid-October. 2.2. Xinanjiang RR model The Xinanjiang RR model was developed by Zhao et al. (1980). It is a lumped RR model (Ju et al., 2009; Li et al., 2009; Zhao, 1992). The model consists of three sub-models, a three layer evapotranspiration sub-model, a runoff generation sub-model and a runoff routing sub-model (Li et al., 2009; Zhao, 1992). The schematic diagram of the model can be found in many references (Li et al., 2009; Zhao, 1992). The model parameters, the state variables, input and output are summarized in Table 2. Issues of complexity and uncertainty of the parameters and input variable uncertainty for the Xinanjinag RR model are very important. Using sensitivity analysis methods, correlation, correlation analysis methods and data assimilation methods to study the state variable and parameter updates are the main topics in this paper. In the Xinanjiang model, EU, EL, ED are linearly correlated with WU, WL, WD, respectively; RS, RSS, RG are linearly correlated with TRS, TRSS, TRG, respectively; R, RB are linearly correlated with FR. So, in this study, the state variables are divided into two categories, the main state variables: WU, WL, WD, FR, S, TRSS, TRG and the secondary state variables: R, RB, W, RS, RSS, RG, TRS, TR. During data assimilation, we only study how to update the main state variables. The other state variables, such as, R, RB, W, RS, RSS, RG, TRS and TR are updated according to the main state variables. 2.3. Particle swarm optimization Particle swarm optimization (PSO) is a population based stochastic optimization technique developed by Kennedy and Eberhart (1995). First the PSO is initialized with a group of random particles (solutions) and then searches for optima by updating each generation. In every iteration, each particle is updated by two ‘‘best’’ values. The first value is the location of the best solution (fitness) a particle has achieved so far. This value is called pBest. The second ‘‘best’’ value is the location of the best solution that any

104

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

Fig. 1. Location map for the study area: (A) The Luo River basin and (B) The Shiguan River basin.

Table 1 Some hydrology and climate statistical characteristics of the two studied catchments.

Luo River basin Shiguan River basin

Area (km2)

Mean of annual rainfall (mm)

Variance of annual rainfall (mm2)

Mean of annual evapotranspiration (mm)

Mean of annual streamflow (m3 s1)

Variance of annual streamflow (m3 s1)2

150 5730

2330 1019

1090 58,103

1606 1000

8.76 352.8

3.406 41,432

Table 2 The parameters, state variables, inputs and outputs of the Xiananjian RR model. (Li et al., 2009). Parameters, variables, input and output Parameters

Evapotranspiration

Runoff generation Runoff routing

State variables

Runoff

Area factor Tension water

Free water Flow

Input Output

Rainfall Pan evaporation Evapotranspiration

Discharge

UM LM DM C B Im Sm Es Kg Ks Cg Ci Cs L R RB RS RSS RG FR W WU WL WD S TRS TRSS TRG TR P Ep E EU EL ED Qsim

Physical meaning

Range and units

Averaged soil moisture storage capacity of the upper layer Averaged soil moisture storage capacity of the lower layer Averaged soil moisture storage capacity of the deep layer Coefficient of the deep layer Exponential of the distribution to tension water capacity Percentage of imprevious and saturated areas in the catchment Areal mean free water capacity of the surface soil layer Exponent of the free water capacity curve influencing the development of the saturated area Outflow coefficients of the free water storage to groundwater relationships Outflow coefficients of the free water storage to interflow relationships Recession constants of the groundwater storage Recession constants of the lower interflow storage Recession constant in the lag and route method Lag in time empirical value Runoff from the previous area Runoff from the imprevious area Surface runoff Interflow runoff Groundwater runoff The runoff producing area proportion Areal mean tension water storage Areal mean tension water storage in the upper layer Areal mean tension water storage in the lower layer Areal mean tension water storage in the deepest layer Areal mean free water storage Surface runoff Interflow Groundwater The total sub-basin inflow to the channel network The measured areal mean rainfall The measured pan evaporation The actual evapotranspirations from the whole basin The evapotranspirations from the upper soil layer The evapotranspirations from the lower soil layer The evapotranspirations from the deepest soil layer The discharge from a sub-basin

5–100 (mm) 50–300 (mm) 5–100 (mm) 0.15 () 0–1 () 0–0.1 (%) 5–100 (mm) 1.2 () 0.05–0.70 () 0.7–0.8 () 0.9–0.999 () 0.05–0.95 () 0.001–0.8 () 0–3(d) m3 s1 m3 s1 m3 s1 m3 s1 m3 s1 – mm mm mm mm mm m3 s1 m3 s1 m3 s1 m3 s1 mm mm mm mm mm mm m3 s1

105

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

neighbor of a particle has achieved so far. This best value is a neighborhood best and is called nBest. When a particle takes the whole population as its neighbors, the best location is a global best and is called gBest, which is the final estimate. PSO performance compares favorably to the performance of various other methods (Kennedy and Spears, 1998), and PSO has been successfully applied for parameter optimization in some hydrology models and environment models (Chau, 2006; Gill et al., 2006). 2.4. Ensemble Kalman filter for state variables calibration The Ensemble Kalman filter (EnKF) is a sophisticated sequential data assimilation method. It applies an ensemble of model states to represent the error statistics of the model estimate, it applies ensemble integrations to predict the error statistics forward in time, and it uses an analysis scheme which operates directly on the ensemble of model states when observations are assimilated. The EnKF has proven to efficiently handle strongly nonlinear dynamics and large state spaces and is now used in realistic applications with primitive equation models for the ocean and atmosphere. Because the EnKF has been described in many papers (Clark et al., 2008; Evensen, 2009; Moradkhani et al., 2005b), the process details are omitted here. 2.5. Ensemble Kalman filter for parameter calibration The hydrological model parameters can be estimated using the priori observation data. However, because of complex error sources, such as forcing data, measurement data, and model structure, it is necessary to re-evaluate the model parameters which should not change with time. A few studies have been performed on simultaneous estimation of the state variables and the model parameters (Moradkhani et al., 2005a; Young, 2002) using filter methods with real-time processing of the state variables in hydrology models and environmental models. In this paper, we first assume that the model parameters do not change with time, and we study how to revise the state variables in the Xinanjiang RR model using EnKF data assimilation. Then we take into account dual data assimilation for the model state variable and the model parameters. A method for jointly estimating the model state variables and parameters is where the state and parameter vectors are concatenated into a single joint state vector (Bras and Restrepo-Posada, 1980; Reichle et al., 2002; Young, 2002). The disadvantage of this approach is that the numbers of the unknown state variables and the parameters and the degrees of freedom are large. This leads to unstable and intractable conditions for non-linear dynamic models. Thus, dual data assimilation becomes the best estimation method for which the state variables and the parameters are estimated simultaneously (Lü et al., 2011; McMillan et al., 2012; Moradkhani et al., 2005a). Dual data assimilation consists of two reciprocal filter processes, the state variable filter process and the parameter filter process. The first dual data assimilation was a dual extended Kalman filter (EKF) (Nelson, 1999; Wan and Nelson, 1997) for estimating neural networks model signals (state) and weights (parameter). The state variables and the parameters are considered separately in the dual EKF method. The two filter processes are parallel. The dual EnKF is a generalization of the dual EKF. The dual EnKF simultaneously takes into account the estimation of the state variables and the parameters. Here we will couple the Xinanjiang RR model and the dual EnKF method to estimate the state variables and the parameters simultaneously and determine the impact of the assimilation time scale on the simulation results. Parameter evolution is set up artificially, i.e., it is assumed that the parameters follow a random walk. Therefore, in EnKF, parameter samples are made as follows:

iþ i hi iþ1 ¼ ht þ st ;

sit  Nð0; Rht Þ

hi tþ1

ð1Þ hiþ t

where is the i th ensemble parameter at time t + 1 and is the i th updated ensemble parameter at time t and sit is noise with covariance Rht . Using the upper parameter sets, a model state ensemble and predictions are made, respectively i iþ i xi tþ1 ¼ f ðxt ; utþ1 ; htþ1 ; tÞ

ð2Þ

i ^itþ1 ¼ hðxi y tþ1 ; hiþ1 Þ

ð3Þ

xi tþ1

where is the ith ensemble member forecast at time t + 1 and xiþ t is the ith updated ensemble member at time t. uitþ1 is input vector at time t + 1, e.g., mean areal precipitation, hi tþ1 is the ith ensemble ^itþ1 is the ith predictive variable at time parameter at time t + 1, and y t + 1, the ‘‘+’’ superscript denotes that the estimate is a posteriori, the ‘‘‘‘ superscript denotes that the estimate is a priori. Updating the parameter ensemble members is made according to the standard Kalman filter equation h i i ^i hiþ tþ1 ¼ htþ1 þ K tþ1 ðytþ1  ytþ1 Þ

ð4Þ

K htþ1

Here, is the Kalman gain for correcting the parameter trajectories and is obtained by 1 yy y K htþ1 ¼ Rhy tþ1 ½Rtþ1 þ Rtþ1 

ð5Þ

hy tþ1

where R is the cross covariance of parameter ensemble and prediction ensemble. 2.6. Sensitivity analysis In this section we investigate the sensitivity analysis of the parameters in the Xinanjiang RR model. Thus, we study how the variation in the output streamflow of the Xinanjiang RR model can be apportioned, qualitatively or quantitatively, to the different model parameters. A sensitivity analysis of the Xinanjiang RR model is important to increase understanding or quantification of the hydrology system (e.g., understanding the relationships between input and output variables); and model development (e.g., searching for errors in the model). Here, the sensitivity analysis is to limit the number of parameters to be updated in the assimilation process. There are several ways, such as local methods, methods based on emulators, screening methods, variance based methods, high dimensional model representations, and methods based on Monte Carlo filtering, to perform sensitivity analyses. In this paper, we use a sampling-based sensitivity method (Helton et al., 2006). By increasing or decreasing input values it is possible to detect the impact of changes on model outputs. Here the changes in model parameters relevant to the Xinanjiang RR model are included. Each parameter is changed from 20% to 20% of their values (denoted: a = 20%, 10%, 10% and 20%, respectively) while the other parameters are fixed. The simulations are carried out for each of the changed values. The following sensitivities are investigated: Daily sensitivity

DSi ¼

Q sim;i ðs þ DsÞ  Q sim;i ðsÞ  100% Q sim;i ðsÞ

ð6Þ

Average sensitivity in the estimation period,

AS ¼

1 XN jDSi j i¼1 N

ð7Þ

Relative sensitivity:

PN 1 jDSi j RS ¼ N Ds i¼1 j s  100j

ð8Þ

106

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

where Qsim,i is simulated streamflow for the ith day, s is input parameter, Ds is a increment of the initial parameter value, and N is the total number of days for the simulation period.

V t þ Rht . According to Liu (2000), the artificial evolution needs to be modified and the conditional evolution density of parameters can be written as follows

2.7. Inputs, outputs and parameter uncertainties

Pðhtþ1 jht Þ  Nðhjtþ1 jahjt þ ð1  aÞht ; h V t Þ

2

The EnKF theory is based on the statistical estimation of sufficient ensemble members. A large ensemble size however requires a large computational resource. There are many methods to determine the ensemble size (Chen et al., 2009; Moradkhani et al., 2005a; Vrugt et al., 2006). It is important to check the plot of the normalized RMSE ratio (NRR) so that an ensemble size can be determined (Moradkhani et al., 2005a; Vrugt et al., 2006). In this study we use an ensemble size of N = 50 following Moradkhani et al. (2005a) and Vrugt et al. (2006), who showed its efficiency in nonlinear hydrological models by studying the normalized RMSE ratio (NRR) of the ensemble size. The spread of the ensemble members is determined by the specified error in the model structure, the input forcing and discharge data. Realistic assumptions for errors in input forcing and responses are essential for proper assimilation of data by filtering. For area average rainfall and evapotranspiration derived in operational flood estimation systems with a limited number of rain gauges, the uncertainties can be up to 50% (Willems, 2001). In this study, a preliminary estimation of input error terms based on Weerts and Serafy (2005, 2006) is used:

Ptrue ¼ Pinput þ dP; dP  Nð0; ð0:15Pinput þ 0:2Þ2 Þ

ð9Þ

and

Etrue ¼ Einput dE

ð10Þ

ð14Þ

where a ¼ 3d1 and d is a factor in (0, 1], which is typically around 2d 0.95–0.99. 2.8. Statistical analysis The root mean square error (RMSE), the Pearson’s correlation (R) and the mean bias error (MBE) of the observed and simulated streamflow are used to assess the performance of the data assimilation scheme. The Nash–Sutcliffe coefficient (Nash and Sutcliffe, 1970) of efficiency (CE) is used to analyze the simulated results. The CE has been widely used to evaluate the performance of hydrologic models. It is the ratio of the mean square error to the variance in the measured data subtracted from unity. It can range from  1 to 1. An efficiency of 1 (CE = 1) corresponds to a perfect match of modeled discharge to the observed data. An efficiency of 0 (CE = 0) indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero (CE < 0) occurs when the observed mean is a better predictor than the model or, in other words, when the model residual (described by the numerator), is larger than the data variance (described by the denominator) (Wilcox et al., 1990). The CE is an improvement over RMSE for model evaluation purposes because it is sensitive to differences in the measured and model estimated means and variances. The Persistence Index (PI) (Kitanidis and Bras, 1980) is used to assess the peak estimate performance of the model.

 X 1=2 N 1 2 ðQ  Q Þ sim;i obs;i i¼1 N

ð15Þ

PN   i¼1 ðQ sim;i  Q sim ÞðQ obs;i  Q obs Þ ffi R ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P PN 2 n 2   i¼1 ðQ sim;i  Q sim Þ i¼1 ðQ obs;i  Q obs Þ

ð16Þ

The daily time series for dE is generated with a uniform random distribution between 0 and 2. In the real flood estimation system this means that the error in evaporation can be 90–100% during the night times. This multiplicative dE ensures that it is possible to account for this effect if this information is in the discharge signal. The discharge produced by the model with input Ptrue and Etrue is called Qtrue. The uncertainty in the discharge measurement can be obtained from the rating curve calibration data for a given gauging station. In this paper, a preliminary estimate of discharge measurement error is chosen similar to that of Georgakakos (1986), who assumed a standard deviation of 0.1 times the measured discharges as given below

RMSE ¼

Q true ¼ Q obs þ dQ; dQ  Nð0; ð0:1Q obs Þ2 Þ

 sim is the averwhere Qsim,i is the simulated streamflow at time i, Q age value of Qsim,i, Qobs,i is the measured streamflow at time i, and  obs is the average of Qobs,i. N is the total number of observations. Q

ð11Þ

The artificial parameter evolution at each time step by adding a small random perturbation provides a new parameter set in the simulations. This method may lose information in a posterior distribution of parameters that are relatively diffuse. West (1993) introduced a kernel smoothing of parameter samples to deal with this problem. Suppose that the h, and their weights wjt , j = 1, . . . , m denote a random measure fhjt ; wjt g to characterize the discrete Monte Carlo approximation to posterior density of parameters as P(ht+1|y1, . . . , yt) with mean  ht and the variance matrix Vt:

Pðhtþ1 jy1 ; . . . ; yt Þ 

m X

2

wjt Nðhtþ1 jmjt ; h V t Þ

ð12Þ

j¼1

where h is the smoothing or variance reduction parameter and

mij ¼ ahjt þ ð1  aÞht ; with a ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1h

ð13Þ

If the Monte Carlo approximation to posterior density P(ht+1|y1, . . . , yt) has mean  ht and variance matrix Vt, the parameter evolution in Eq. (13) with independent perturbation hjt , which was assumed to be independent of ht, has the correct mean  ht but variance matrix

MBE ¼

1 XN ðQ sim;i  Q obs;i Þ i¼1 N PN

CE ¼ 1  Pi¼1 N

ðQ obs;i  Q sim;i Þ2  obs Þ2 ðQ Q

i¼1

ð17Þ

ð18Þ

obs;i

3. Results and discussion 3.1. Parameter optimization The PSO algorithm is linked to the Xinanjiang RR model in a MATLAB environment so that the PSO minimizes the error function calculated by comparing observed and simulated hydrographs. Thirteen model parameters (all except L) are involved in the PSO process as independent variables. The PSO algorithm with fifty particles and one hundred epochs was performed in such a way that each independent variable could vary in its feasible domain (variable range). Rainfall-runoff events were recorded from 1988 to 2003 in the Luo river basin, and flood hydro-graphs were recorded at the Nangao hydrometric station. Thirteen model parameters were considered as independent variables to minimize the calibration objective function which is the sum of absolute errors using rainfall-runoff data from 1988 to 1999. Rainfall data represent an

107

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114 Table 3 The parameter analysis of the Xinanjiang RR Model, here the median values are initial values for Cases 5–6. Parameter

Sensitivity

Optimization value

Median

UM LM DM B Im C Sm Es Kg Ki Cg Ci Cs L

No No No Yes No No Yes No Yes No No Yes Yes –

20 60 40 0.3 0.02 0.15 40 1.2 0.39 0.31 0.992 0.7 0.3 0

52.5 175 52.5 0.5 0.05 0.075 52.5 1.2 0.375 0.375 0.9495 0.5 0.4005 0

3.2. State variables

average of four raingauges. Table 3 shows the optimized results of the model parameters using the 1988–1999 data. The calibration error (RMSE) is 31.0 m3 s1. It is not credible to use the model parameters optimized with historical data to predict the streamflow into the future because the mean of annual stramflow is only 8.76 m3 s1. The model parameters are uncertain and cannot ensure the precision and reliability of the predicted results. In reality,

In the Xinanjiang model, all of the state variables are shown in Table 2. The main state variables are WU, WL, WD, FR, S, TRSS and TRG. Because of the strong nonlinearity of the state-to-output transition equations, the EnKF linear updating rule can only be used to update a subset of the main state variables. So we must evaluate the correlation between the main state variables and the output streamflow and select a subset of the main state variables to update with the EnKF data assimilation. The other state variables can change, respectively. On the other hand, Qsim(t) has a strong relation with Qsim(t  1). So we need to consider the correlation relation between Qsim(t  1) and Qsim(t). In the correlation relation analysis process, we model the output streamflow of 2000–2003 using the following condition: the model parameters are optimization values (Table 3) and the initial state variables are: WU = 20; WL = 50; WD = 10; FR = 0.89; S = 2; TRSS = 0.3Qobs; TRG = 0.4Qobs Qsim|t=0 = Qobs. The Pearson’s correlation (R) between the state variables and the simulated output streamflow was calculated with four years of simulated data. Fig. 2 shows two-dimensional scatter

1

1

0.8

0.8

R=0.4527

WLnorm[-]

WUnorm[-]

the model parameters and state variables are random values with varying degrees of uncertainty. These uncertainties may have a profound effect on the simulated results.

0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

R=0.1826

FRnorm [-]

WDnorm[-]

0.8 0.6 0.4 0.2 0

0.2

0.4

0.6

0.8

TRSSnorm[-]

Snorm[-]

R=0.8226

0.6 0.4 0.2 0.2 +++++ +

0.4

0.6

0.8

0.4

0.6

0.8

1

R=0.5204

0.6 0.4 0.2 0

0.2

0.4

0.6

0.8

1

0.8 0.6 0.4 0.2 0.2

R=0.8433

0.6 0.4 0.2 0

0.2

0.4

0.6

0.8

1

1

++

R=0.5876

0

0.8

0

1

Qsim(t-1)norm [-]

0

1

TRGnorm[-]

0.2

1

0.8

0

0

0.8

0

1

1

0

0.2

1

1

0

0.4

0

1

R=0.3883

0.6

0.4

0.6

Qsimnorm [-]

0.8

1

0.8

R=0.7868

0.6 0.4 0.2 0

0

0.2

0.4

0.6

0.8

1

Qsimnorm [-]

Fig. 2. Two-dimensional scatter plots of relationships between the normalized state variables and the output streamflow at Luo River basin (Year 2000–2003).

108

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

plots which describe the correlation relations and correlation coefficients of the output streamflow Qsim at time t and the state variables at time t for four years of simulation results (year 2000–2003). The horizontal axes are normalized streamflow values at time t and the vertical axes are the normalized state variables at time t. The figure indicates that some of the state variables have strong linear correlations with the output streamflow. The four largest correlation coefficients were obtained for the TRSS, S, Qsim (t  1) and TRG. Therefore, we decided to recursively update only the following model variables: S, TRSS, TRG and Qsim (t  1). On the other hand, many values of the WU, WL and WD equal to 1. The main reason is that the Xinanjiang model is based on the stored-full runoff theory, and WU, WL, WD reflect areal mean tension water storage in the upper layer, the lower layer and the deepest layer. When rain exceeds a certain value, WU, WL, WD reached their maximum values, and WU, WL and WD do not change further as rainfall increases. FR is the runoff producing area proportion. The reason that FR values do not exceed 0.6 is that the runoff producing area proportion is 60% during the simulation period.

3.3. Parameter sensitivity Parameter ranges in the Xinanjiang RR model are shown in Table 2. Some of the parameters are sensitive, and small changes may cause large changes in the simulated results. Other parameters have little impact on the model results. Among the Xinanjiang RR model parameters, L is an integer, equal to 0, 1, 2 or 3 (Zhao, 1992). According to the characteristics of the study area (the basin

Relative sensitivity [-]

3.0x10

3.4. Data assimilation 3.4.1. Presentation of the study cases In this paper, we study seven cases (Table 4). In Cases 1–4, we assume that the model parameters can be obtained by optimization with the historical data (1988–1999). We use these results to simulate the rainfall-runoff processes of 2000–2001. In the simulated processes, we do not include data assimilation (Case 1); do include parameter assimilation (Case 2); only include state variable assimilation (Case 3) and include simultaneously parameter assimilation and state variable assimilation (Case 4). To apply the proposed model in an ungauged basin (it is not exactly an

-03

-03

3.0x10

Um Lm Dm C

B Im

-03

2.5x10

2.5x10

-03

2.0x10

-03

2.0x10

1.5x10

-03

1.5x10

1.0x10

-03

1.0x10

5.0x10

-04

5.0x10

0.0x10

+00

0.0x10

-03

-03

-03

-04

-0.2

+00

-0.1

0

0.1

-03

2.5x10

-0.2

0.2

Sm Kg Ki Eg

-03

3.0x10

Relative sensitivity [-]

area is relatively small), we set L = 0 (constant) in this study. The sensitivity of the other parameters was studied. For the four different changes (20%, 10%, 10%, and 20%) of the initial values of each the Xinanjiang RR model optimization parameters, 13 curves were obtained by carrying out the simulations using four years of forcing data (2000–2003) (Fig. 3). The relative sensitivities were simultaneously computed. Then sums of the relative sensitivities and the order of the parameters from large to small: Cs, Ci, Kg, B, Sm, Ki, Dm, Um, Lm, Cg, Eg, Im and C. If the sum is larger than 0.0079 (for Sm), we identify the parameter as sensitive. Most of the parameters are not sensitive. Table 3 shows the relative sensitivities of the 13 parameters. According to the results, the parameters: Cs, Ci, Kg, B and Sm are adjusted in the parameter data assimilation process. The other parameters do not change with time, and their values used in the simulations are presented in Table 3.

-0.1

0

-01

0.1

0.2

Ci Cs Cg

1.0x10

-02

8.0x10

-03

2.0x10

-02

6.0x10

-03

1.5x10

-02

4.0x10

-03

1.0x10

-02

2.0x10

-04

5.0x10

+00

0.0x10

0.0x10 -0.2

-0.1

0

α [-]

0.1

0.2

+00

-0.2

-0.1

0

0.1

0.2

α [-]

Fig. 3. Relative sensitivity of the Xinanjiang RR model parameters at Luo River basin (Year 2000–2003). The Y-scale for the first three graphs has the same range, but the Y-scale on the fourth graph represents differentials.

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114 Table 4 Some cases studied in this paper. OP: Optimization; PA: Parameter assimilation; SA: State variable assimilation; TS: Time scale for assimilation; ME: Median. Median values are shown in Table 3. TS = 3 days for all cases. Case Luo River basin

Shiguan River basin

1 2 3 4 5 6 7 1 3 4

Initial parameter

PA

SA

OP OP OP OP ME ME ME OP OP OP

No Yes No Yes Yes No Yes No No Yes

No No Yes Yes No Yes Yes No Yes Yes

ungauged basin), we only assume that no historical data are available and that the model parameters cannot be determined using the proposed PSO algorithm. In this case, the model parameters are set at the median values (Table 3) and the initial state variables are stochastic values selected from the ranges of values. The median values are computed by averaging the sum of the minimum and maximum values of the parameters. In addition, for only state variable assimilation, only parameters assimilation, or dual data assimilation, we study another three cases (Cases 5–7). 3.4.2. Results without data assimilation In Case 1, the model parameters are estimated using the PSO algorithm, then the Xinanjiang RR model is run without any data assimilation. The simulated results (streamflow, Qsim) are shown in Fig. 4. For the sake of analysis, the whole rainfall time series in the simulated period (2000–2001) is also shown in Fig. 4. The total rainfall of the entire simulation period (Year: 2000–2001) is 5968 mm. In Case 1, the simulated results were not satisfactory. In the overall analysis period, RMSE was large (RMSE = 15.36 m3 s1) and the correlation coefficient (R) was relatively small (R = 0.868). The simulated values overestimated the observed values. Especially, at DOY170, DOY199 and DOY 244, the percents for overestimated values are 94.8%. 77.8% and 67.6%. The reasons are (1) the observed precipitation data are large; (2) the parameter errors are large. For accurate flood estimation, it is very important to accurately estimate the streamflow during wet periods. To analyze the simulated results thoroughly, we divided the observed streamflow into three periods (Table 6), dry period (December, January and February), level period (March, April, May, October and November) and wet period (May, June, July, August and September). In each period, the RMSE was consistently large with a maximum value of 23.79 m3 s1 for the wet period and a minimum value of 2.12 m3 s1 for the level period. For the wet period, the simulated values completely overestimated the streamflow. During the level period, the simulation results were favorable with correlation coefficients reaching 0.906. 3.4.3. Results in a gauged basin with data assimilation Fig. 4 shows the simulated results Case 2 (OP + PA), Case 3 (OP + SA) and Case 4 (OP + PA + SA) (see Table 4, too). In the discharge for Case 2 (OP + PA), the measured value of the flood peak is larger than the simulated value, and the simulated result is obviously later than the measured recession time. However, because of the parameter calibration, the simulation results during DOY150 and DOY250 are better than that in Case 1 (OP). So, parameter assimilation improves the simulation of the flood peak. There is a problem in the streamflow chart of Case 2 in that the simulation slightly overestimates the streamflow for the recession flow. The measured values in Case 3 are smaller than the simulated values

109

for the recession flow. Although some of the recession curves do match well, the simulated recession curve is a little lower than the measured curve. For DOY 200, the simulation is similar to Case 1, in that the simulation value seriously overestimates the streamflow. However, the simulation in Case 2 is better than that in Case 1. From the global aspect, the simulated discharge curves of Cases 2 and 3 are better than that of Case 1. Observing the hydrograph of Case 4, it avoids overestimating the flood peak during DOY 150 and DOY 250 and the recession flow. For Case 3 (OP + SA) in Fig. 4, the residual error is large when the discharge is large, and the changes in the discharge and residual curves are somewhat synchronized. The values of the residuals are stable when the discharge is small, and most of them are negative except when there is a sudden change at the flow peak. For Case 4 (OP + PA + SA), the residual values fluctuate, but the range change is much smaller than for Cases 2 (OP + PA) and 3 (OP + SA). Throughout the simulation period, the RMSE for Case 4 is only 1.06 m3 s1. For Cases 2 and 3, RMSE values are greater than 10 m3 s1. The RMSEs are less than in Case 1. For Case 4 the correlation coefficient reached 0.999. These results confirm that dual data assimilation is effective. The effectiveness of dual data assimilation is apparent for the various streamflow periods (Table 6). For the wet period, the predicted results are very accurate. For Cases 2, 3 and 4, the RMSE’s are 16.1, 16.1, and 1.2 m3 s1. The correlation coefficients are 0.787, 0.900 and 0.999. For the level period, the RMSE’s are 3.5, 1.9 and 1.1 m3 s1. It can be observed for the error statistical analysis (Table 5) that the parameter alone data assimilation (Case 2) is less than the state variable alone data assimilation (Case 3). During the dry period and the level period, the differences between the assimilation methods are not great. For Case 4 (OP + PA + SA), the simulation results are accurate during every simulated period. 3.4.4. Results in an ungauged basin with data assimilation Fig. 4 shows the simulated discharge results of Cases 5 (ME + PA), 6 (ME + SA) and 7 (ME + PA + SA) too. In these three cases, the selected model is applied in the Luo River basin that was considered as ungauged for parameter calibration. We assume that there is not enough historical data to optimize the model parameter values, so the model parameters are assigned as median values. The initial state variables are similar to those in Case 1. In Case 5 (ME + PA), the Xinanjiang RR model was coupled with the parameter assimilation. In Case 6 (ME + SA), the Xinanjiang RR model was coupled with the state variable data assimilation. In Case 7 (ME + PA + SA), the Xinanjiang model was fused with the dual data assimilation. It can be seen from the Case 5 (ME + PA) simulated hydrograph that when the parameters are median values, the simulation overestimated the observed streamflow during the first year (DOY1365). The reasons may be due to the effects of parameters which overestimate the production flow area or underestimate the evaporation. In the second year (DOY365-730), the simulation is improved. Especially for DOY550-DOY650, the simulation does not overestimate the streamflow. It is obvious that the parameters are calibrated during the parameter data assimilation. For Case 6 (ME + SA), the simulated values are similar to the observed values except for DOY 200. The simulated value is extremely high on DOY 200. The reason for the extreme events at DOY 200 is likely that the rainfall data are inaccurate. For Case 7 (ME + PA + SA), parameter assimilation and state variables assimilation provides accurate simulations for the whole period. The simulation value is similar to the measured value on DOY 200. From the statistical error analysis, we can see that simulation results of Case 5 (ME + PA) are worse than for the other cases (Case 6 and 7). For Case 5 the RMSE is 13.39 m3 s1, CE is only 0.473. The simulation is better than for Case 1 (OP). So, if the initial parameters

Streamflow [m3/s]

300

Streamflow [m3/s]

300

Streamflow [m3/s]

300

Streamflow [m3/s]

300

Streamflow [m3/s]

300

Streamflow [m3/s]

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

300

Streamflow [m3/s]

110

300

Observation Simulation

Case 1 : OP

200 100 0

0

100

200

300

400

500

600

700

800

200

300

400

500

600

700

800

200

300

400

500

600

700

800

200

300

400

500

600

700

800

200

300

400

500

600

700

800

200

300

400

500

600

700

800

300

400

500

600

700

800

Case 2 : OP+PA 200 100 0

0

100

Case 3 : OP+SA

200 100 0

0

100

Case 4 : OP+PA+SA

200 100 0

0

100

Case 5 : ME+PA

200 100 00

100

Case 6 : ME+SA

200 100 0

0

100

Case 7 : ME+PA+SA

200 100 0

0

100

200

Day of Year (2000-2001) Fig. 4. The streamflow simulations for Case 1–7 for 2000–2001.

are medians, and the parameters are calibrated during the simulation period using data assimilation algorithm, the simulation results are better than for the case using historical data to calibrate

parameters with the parameters remaining unchanged during the simulation period (Case 1). The simulation results of Case 6 (ME + SA) are better than those of Case 5 (ME + PA). This result

111

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114 Table 5 Statistical characteristics of the Luo River basin (the flood period: 2000–2001) and the Shiguan River basin (the flood period: 1999). R ()

MBE (m3 s1)

Luo River basin

CE ()

1 2 3 4 5 6 7

15.36 10.56 10.44 1.06 13.39 9.51 4.80

0.8681 0.8199 0.9118 0.9988 0.9063 0.9031 0.9764

0.85 0.04 0.69 0.52 0.63 0.53 0.38

0.3062 0.6721 0.6792 0.9967 0.4727 0.7338 0.9323

Shiguan River basin

1 3 4

31.36 27.80 26.99

0.7827 0.8581 0.8852

8.58 8.15 9.68

0.5722 0.6639 0.6831

350 300

Simulation [m3/s]

RMSE (m3 s1)

Case

OP OP+SA OP+SA+PA

400

mirrors the pre-group (Case 2, 4) result. It shows that the impact of the state variable data assimilation process is larger than that of the parameter data assimilation process. On the other hand, the results in the ungauged case with SA (Case 6) are better than those in the gauged case with SA (Case 3). This is probably obtained by chance here and the reverse could be expected in the general case. In Case 7 (ME + PA + SE), because of the simultaneous data assimilation, the simulation is very accurate (Table 5). The statistical error analysis values are similar to Case 4 (OP + PA + SA) values. The Case 7 (ME + SA + PA) results are superior to those of Cases 2–3 and Cases 5–6. When the parameters are set to their median values, the simulated results of the Xinanjiang RR model coupled with dual data assimilation are relatively accurate. The EnKF of the parameters and the state variables has a very important influence on the real-time calibration of the model parameters. When the model parameters are median and are not real-time calibrated using EnKF data assimilation, only the state variables are real-time calibrated, the errors are relatively large (Case 6). Meanwhile, for the parameter only calibration using the EnKF data assimilation, the simulated results are relatively poor. Table 6 shows the statistical error analysis results for three cases (Case 5–7) during three simulation periods, dry period, level period and wet period. In the dry period, because the runoff is small, the error for the three cases is similar. The RMSE is between 1.73 and 2.57 and R is between 0.821 and 0.888. In the level period, it is obvious that Case 6 (ME + SA) and Case 7 (ME + PA + SA) are better than Case 5 (ME + PA). In the wet period, Case 7 (ME + PA + SA) is the best. The RMSE is only 6.97 m3 s1 and the RMSE is 20.74 m3 s1 for Case 5 (ME + PA) and 14.64 m3 s1 for Case 6 (ME + SA). 3.4.5. Results in Shiguan River basin with data assimilation Fig. 5 shows the simulated streamflow of the Shiguan river basin using the proposed method, and Table 5 shows the statistical analysis for the simulated streamflow of the Shiguan River basin. Here we report simulation results for three kinds of cases (OP,

250 200 150 100 50 0

0

50

100

150

200

250

300

350

400

Observation [m3/s] Fig. 5. The streamflow simulation for Shiguan River basin.

OP + SA, OP + SA + PA) in the Shiguan river basin. The results show that the simulations using the OP + SA and the OP + SA + PA methods are better than the results of OP. RMSE is 31.36 m3 s1 for the OP method, 27.80 m3 s1 for the OP + SA method and 26.99 m3 s1 for the OP + SA + PA method. The Pearson’s correlation is 0.7827 vs. 0.8581 vs. 0.8852 for OP, OP + SA and OP + SA + PA. The OP + SA + PA method is slightly better than the OP + SA. Especially, for the simulation of flood peak in DOY 179, the measured peak for 380 m3 s1, the simulation value for 134 m3 s1 for OP, 273 m3 s1 for OP + SA and 310 m3 s1 for OP + SA + PA. For the large river basin, the streamflow simulation results also have improved by using the proposed method.

3.5. Evolution of model parameters In this section we analyze the evolution of the Xinanjiang RR model parameters in three cases of data assimilation. Many hydrologists have pointed out that estimates of hydrologic model parameters are generally error-prone because past data (e.g., rainfall, streamflow) are used to optimize the model. Real-time calibration of the model parameters is very important. In this paper, the model parameters are adjusted by an automatic algorithm (EnKF algorithm) between each run of the model and the corresponding observation until an optimal parameter set has been found. Fig. 6 shows the optimization results and the evolution of five model

Table 6 The error statistic is computed for different flow groups at the Luo River basin flood examination period (2000–2001). Dry period: December, January, February; Level period: March, April, May, October, November; Wet period: May, June, July, August, September. Flow group

Case Case Case Case Case Case Case

1 2 3 4 5 6 7

Dry period

Level period

Wet period

RMSE (m3 s1)

R ()

RMSE (m3 s1)

R ()

RMSE (m3 s1)

R ()

2.67 1.75 1.78 0.71 1.73 2.51 2.57

0.8963 0.6758 0.871 0.9729 0.8881 0.8211 0.8252

2.12 3.54 1.93 1.14 2.31 1.54 2.02

0.9065 0.6748 0.9388 0.9815 0.8951 0.9576 0.9274

23.79 16.10 16.14 1.17 20.74 14.64 6.97

0.856 0.7869 0.8998 0.9992 0.901 0.887 0.9759

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

Sm (mm)

112

50 40 30

100

Kg (-)

1

B (-)

300

400

500

600

700

800

300

400

500

600

700

800

300

400

500

600

700

800

300

400

500

600

700

800

300

400

500

600

700

800

Parameter Kg

100

1

200

Parameter B

0.5

0

100

1

Ci (-)

200

0.5

0

200

Parameter Ci

0.5

0

100

1

Cs (-)

Case1:OP Case2:OP+PA Case5:ME+PA

Parameter Sm

200

Parameter Cs

0.5 0

100

200

Day of Year (2000-2001) Fig. 6. The evolution of parameters in the data assimilation process for Cases 2 and 5.

parameters in Case 2 (OP + PA) and Case 5 (ME + PA). Because the parameter evolution in Case 4 (OP + PA + SA) and Case 7 (ME + PA + SA) is similar to that in Case 2 (OP + PA) and Case 5 (ME + PA), we do not present detailed results. In Case 2 (OP + PA), the initial values of the model parameters are the optimized values. In the EnKF data assimilation process, the change of the parameters is large during the initial stage. Parameter Kg and B had relatively large changes compared to the other parameters. As simulation time progresses, all parameters gradually approach constant values. For Case 5 (ME + PA), because the initial values deviated from the optimal values during the initial 100 days, several parameters changed slowly (see error bar in Fig. 6). The parameters approached the optimal values or fixed values. This indicates that some optimal results are reliable and some optimal results are not the best during the estimation process. The model parameter data assimilation process, even with the initial error, can correct the model parameters over time. 3.6. Impact of a lag in availability of input data In this section, a combination of the Xinanjiang RR model and the EnKF data assimilation is devised in such a way that updates are not done throughout the fixed time series. We assume that

the precipitation and evapotranspiration is available daily, but the streamflow data are not available at all times. The assumption is that the streamflow data become available in n-day sets (it is not available for the first n  1 days, it is available at the next day, the nth day, and so on). This is a very realistic assumption and can be used to see the improvements in the model performance under conditions of a lag in availability of input data. In this paper, we suppose that n = 2,3,4 or 5. We study the responses of the simulation results on the data assimilation time scale in Cases 2, 3 and 4. When the observations become available, the state variables or parameters are updated. These updated predictions replace the streamflow value in the inputs. The other conditions for Cases 2– 4 do not change. It is found that there are no significant impacts on RMSE and CE. Only for Case 3 (OP + SA), RMSE has an increase when the assimilation time scale changes from 2 to 5 days. For Cases 2 (OP + PA) and 4 (OP + PA + SA), RMSE is almost a constant. From the analysis of CE, Case 4 (OP + PA + SA) is an accurate and robust model. The CE values are close to 1 for different assimilation time scales (n = 2, 3, 4). In Cases 2 (OP + PA) and 3 (OP + SA), the CE changes as the data assimilation time scale changes. It is very sensitive to assimilation time scale in Case 3 (OP + SA). These main results indicate that the assimilation time scales did not contribute to a significantly improved model for the spatial distribution of

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

daily runoff in Case 4 (OP + PA + SA). Case 4 (OP + PA + SA) is a robust model for making daily runoff estimation.

4. Conclusions Although the lumped RR model is widely applied to predict streamflow, due to nonlinearities and frequent human activities, uncertainty in the model parameters forces data uncertainty, and the traditional method to determine model parameters using historical data is constrained. Global intelligence stochastic optimization methods can help to determine model parameters. Recent developments of ensemble-based techniques and filtering strategies for sequential data assimilation now make it possible to implement techniques that properly account for data (input and output), state variables, parameters, and model structural uncertainties, although the details of how such implementations can be achieved are still an active area of research. In this article, the Xinanjiang RR model is applied to estimate streamflow in a southern region of China. The following problems are addressed: the sensitivity of the Xinanjiang RR model parameters, correlation between the state variables and the model output streamflow, determining the model parameters using the PSO algorithm; the real time calibration of state variables using EnKF data assimilation, the real time calibration of the model parameters using EnKF data assimilation, the dual data assimilation for the state variables and the model parameters, the evolution of the model parameters under the data assimilation and the impact of the time step on the data assimilation. Using the PSO algorithm to determine the model parameters and two years of the rainfall-runoff data as testing data, due to model uncertainty, the Xinanjiang RR model are combined with the EnKF method to adjust the state variables or the model parameters or both simultaneously in a real-time estimation process. The results show that coupling the PSO and the dual EnKF data assimilation is very effective for calibrating the simulated streamflow. When the model parameters are optimized with the PSO algorithm, the simulated results are not significantly improved in Case 3 (OP + SA). The simulated results are significantly improved in Case 2 (OP + PA) in which the model parameters are assimilated. The Xinanjiang RR model combined with the EnKF method is applied to an ungauged basin. We assume that the Xinanjiang RR model parameters cannot be determined or that the historical data are scarce, and the initial model parameter values are set as their median values. We compared the simulated results for three cases (the Xinanjiang model coupled with the parameter assimilation or the state variable assimilation, or the dual data assimilation). The dual data assimilation results were better than those for the other cases. The parameter only assimilation was better than the state variable only assimilation results. The model parameters were more important than the state variable assimilation. We also studied the responses of simulated results to different assimilation time scales, 2 days, 3 days, 4 days and 5 days in Cases 2 (OP + PA), 3 (OP + SA) and 4 (OP + PA + SA). The simulated results were very sensitive to the assimilation time scale in Case 3 (OP + SA). The simulated results were not sensitive to the assimilation time scale in Case 4 (OP + PA + SA). Case 4 (OP + PA + SA) is an accurate and robust model for making daily runoff estimation. Through this study, we can see, for the lumped hydrological model (such as Xinanjiang RR model), some internal state variables are strongly correlated with the runoff, (for example TRSS, S and TRG) and some parameters are sensitive (for example Cs, Ci, Kg, B, Sm). When there is observed streamflow data, the data assimilation using the proposed method for the strongly correlated state variables and the strongly sensitive parameters can appropriately improve the simulation results and advance the simulation accuracy.

113

Acknowledgements This research is supported by the NNSF of China (51190090) and the National Basic Research Program of China (2010CB951101) and the NNSF of China (50939006), the special funding of the key laboratory (1069-50985512).

References Akhtar, M., Ahmad, N., Booij, M.J., 2008. The impact of climate change on the water resources of Hindukush–Karakorum–Himalaya region under different glacier coverage scenarios. J. Hydrol. 355, 148–163. Bingeman, A., Kouwen, N., Soulis, E.D., 2006. Validation of the hydrological processes in a hydrological model. J. Hydrol. Eng. 11, 451–463. Bras, R., Restrepo-Posada, P., 1980. Real time automatic parameter calibration in conceptual runoff forecasting models. In: Proceedings of the Third International Symposium on Stochastic Hydraulics, pp. 61–70. Burnash, R.J.C., Ferral, R.L., McGuire, R.A., 1973. A generalized streamflow simulation system-conceptual modeling for digital computers. National weather service and California department of water resources. Chau, K.W., 2006. Particle swarm optimization training algorithm for ANNs in stage prediction of Shing Mun River. J. Hydrol. 329, 363–367. Chen, C., Malanotte-Rizzoli, P., Wei, J., Beardsley, R.C., Lai, Z., Xue, P., Lyu, S., Xu, Q., Qi, J., Cowles, G.W., 2009. Application and comparison of Kalman filters for coastal ocean problems: an experiment with FVCOM. J. Geophys. Res. 114, C05011. Cheng, C.T., Ou, C.P., Chau, K.W., 2002. Combining a fuzzy optimal model with a genetic algorithm to solve multi-objective rainfall-runoff model calibration. J. Hydrol. 268, 72–86. Clark, M.P., Rupp, D.E., Woods, R.A., Zheng, X., Ibbitt, R.P., Slater, A.G., Schmidt, J., Uddstrom, M.J., 2008. Hydrological data assimilation with the ensemble Kalman filter: use of streamflow observations to update states in a distributed hydrological model. Adv. Water Resour. 31, 1309–1324. Duan, Q., Gupta, V.K., Sorooshian, S., 1993. A shuffled complex evolution approach for effective and efficient global minimization. J. Optim. Theory Appl. 76, 501– 521. Duan, Q., Sorooshian, S., Gupta, V.K., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28, 1015–1031. Evensen, D., 2009. Data assimilation-the ensemble Kalman filter, second ed. Springer. Evensen, G., 1994. Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 99, 10143–10162. Gan, T.Y., Dlamini, E.M., Biftu, G.F., 1997. Effects of model complexity and structure, data quality, and objective functions on hydrologic modeling. J. Hydrol. 192, 81–103. Georgakakos, K.P., 1986. A generalized stochastic hydrometeorological model for flood and flash flood forecasting: 2. Case studies. Water Resour. Res. 22, 2096– 2106. Gill, M.K., Kaheil, Y.H., Khalil, A., McKee, M., Bastidas, L., 2006. Multiobjective particle swarm optimization for parameter estimation in hydrology. Water Resour. Res. 42, W07417. Helton, J.C., Johnson, J.D., Salaberry, C.J., Storlie, C.B., 2006. Survey of sampling based methods for uncertainty and sensitivity analysis. Reliab. Eng. Syst. Saf. 91, 1175–1209. Jayawardena, A.W., Zhou, M.C., 2000. A modified spatial soil moisture storage capacity distribution curve for the Xinanjiang model. J. Hydrol. 227, 93–113. Jones, J.P., Sudicky, E.A., McLaren, R.G., 2008. Application of a fully-integrated surface-subsurface flow model at the watershed-scale: a case study. Water Resour. Res. 44 (March), W03407. Ju, Q., Yu, Z., Hao, Z., Ou, G., Zhao, J., Liu, D., 2009. Division-based rainfall-runoff simulations with BP neural networks and Xinanjiang model. Neurocomputing 72, 2873–2883. Kalman, R., 1960. A new approach to linear filtering and prediction problems. Trans. ASME J. Basic Eng. 82, 35–45. Kennedy, J., Eberhart, R.C., 1995. Particle swarm optimization. In: Proceedings of IEEE International Conference on Neural Networks, vol. IV. pp. 1942–1948. Kennedy, J., Spears, W.M., 1998. Matching algorithms to problems: An experimental test of the particle swarm and some genetic algorithms on the multimodal problem generator. In: Proceedings of the IEEE Congress on Evolutionary Computation. pp. 78–83. Kitanidis, P.K., Bras, R.L., 1980. Real-time forecasting with a conceptual hydrologic model. Part II: applications and results. Water Resour. Res. 16, 1034–1044. Kuczera, G., 1998. On validity of first-order prediction limits for conceptual hydrological models. J. Hydrol. 103, 229–247. Li, H., Zhang, Y., Chiew, F.H.S., Xu, S., 2009. Predicting runoff in ungauged catchments by using Xinanjiang model with MODIS leaf area index. J. Hydrol. 370, 155–162. Li, X., Contreras, S., Sole-Benet, A., 2007. Spatial distribution of rock fragments in Dolines: a case study in a semiarid Mediterranean mountain-range. Catena 70, 366–374.

114

H. Lü et al. / Journal of Hydrology 480 (2013) 102–114

Lindstrom, G., Johannson, B., Persson, M., Gardelin, M., Bergstrom, S., 1997. Development and test of the distributed HBV-96 hydrological model. J. Hydrol. 201, 272–288. Liu, D., Yu, Z., Hao, Z., Yang, C., Ju, Q., 2007. Groundwater simulation in the Yangtze River basin with a coupled climate-hydrologic model. J. China Univ. Geosci. 18, 155–157. Liu, F., 2000. Bayesian time series: analysis methods using simulation based computations. Master’s Thesis. PhD Dissertation. Institutes of Statistics and Decision Sciences, Duke University. Liu, Y.Q., Gupta, H.V., 2007. Uncertainty in hydrologic modeling: toward an integrated data assimilation framework. Water Resour. Res. 43, W07401. Lü, H., Li, X., Yu, Z., Horton, R., Zhu, Y., Hao, Z., Xiang, L., 2010. Using a h-1 filter assimilation procedure to estimate root zone soil water content. Hydrol. Process. 24, 3648–3660. Lü, H., Yu, Z., Zhu, Y., Drake, S., Hao, Z., Sudicky, E.A., 2011. Dual state-parameter estimation of root zone soil moisture by optimal parameter estimation and extended Kalman filter data assimilation. Adv. Water Resour. 34, 395–406. Lü, H., Zhu, Y., Skaggs, T.H., Yu, Z., 2009. Comparison of measured and simulated water storage in dryland terraces of the Loess Plateau, China. Agric. Water Manag. 96, 299–306. McCabe, M.F., Franks, S.W., Kalma, J.D., 2005. Calibration of a land surface model using multiple data sets. J. Hydrol. 302, 209–222. McMillan, H.K., Hreinsson, E.Ö., Clark, M.P., Singh, S.K., Zammit, C., Uddstrom, M.J., 2012. Operational hydrological data assimilation with the retrospective ensemble Kalman filter: use of observed discharge to update past and present model states for flow forecasts. Hydrol. Earth Syst. Sci. Discuss. 9, 9533–9575. Misirli, F., Gupta, H.V., Sorooshian, S., Thiemann, M., 2003. Bayesian recursive estimation of parameter and output uncertainty for watershed models: calibration of watershed models. Water Sci. Appl. Ser. 6, 113–124. Moore, R.J., 2007. The PDM rainfall-runoff model. Hydrol. Earth Syst. Sci. 11, 483– 499. Moradkhani, H., Hsu, K.-L., Gupta, H., Sorooshian, S., 2005a. Uncertainty assessment of hydrologic model states and parameters: sequential data assimilation using the particle filter. Water Resour. Res. 41, W05012. Moradkhani, H., Sorooshian, S., Gupta, H., Houser, P.R., 2005b. Dual state-parameter estimation of hydrological models using ensemble Kalman filter. Adv. Water Resour. 28, 135–147. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models. Part I – a discussion of principle. J. Hydrol. 10, 282–290. Nelson, A.T., 1999. Nonlinear estimation and modeling of noisy time series by dual Kalman filtering methods. Ph.D. Thesis. Oregon Graduate Institute. Reichle, R.H., Mclaughlin, D.B., Entekhabi, D., 2002. Hydrologic data assimilation with the ensemble Kalman filter. Mon. Weather Rev. 130, 103–114.

Salamon, P., Feyen, L., 2009. Assessing parameter, precipitation, and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter. J. Hydrol. 376, 428–442. Serbert, J., 2000. Multi-criteria calibration of a conceptual runoff model using a genetic algorithm. Hydrol. Earth Syst. Sci. 4, 215–224. Thiemann, M., Trosset, M., Gupta, H., Sorooshian, S., 2001. Bayesian recursive parameter estimation for hydrologic models. Water Resour. Res. 37, 2521– 2535. Vrugt, J.A., Gupta, H.V., Nualláin, B., Bouten, W., 2006. Real-time data assimilation for operational ensemble streamflow forecasting. J. Hydrometeorol. 7, 548–564. Wan, E.A., Nelson, A.T., 1997. Dual Kalman filtering methods for nonlinear prediction, smoothing, and estimation. In: Advances in Neural Information Processing, vol. 9. Weerts, A.H., Serafy, G.Y.H.E., 2005. Particle filtering and ensemble Kalman filtering for input correction in rainfall runoff modeling. In: International Conference on Innovation Advances and Implementation of Flood Forecasting Technology, 17– 19 October 2005. Troms, Norway, 1–15. Weerts, A.H., Serafy, G.Y.H.E., 2006. Particle filtering and ensemble Kalman filtering for state updating with hydrological conceptual rainfall-runoff models. Water Resour. Res. 42, W09403. West, M., 1993. Mixture models, Monte Carlo, Bayesian updating and dynamic models. Comput. Sci. Statist. 24, 325–333. Wilcox, B.P., Rawls, W.J., Brakensiek, D.L., Wright, J.R., 1990. Predicting runoff from rangeland catchments: a comparison of two models. Water Resour. Res. 26, 2401–2410. Willems, P., 2001. Stochastic description of the rainfall input errors in lumped hydrological models. Stoch. Environ. Res. Risk Assess. 15, 132–152. Xie, X., Zhang, D., 2010. Data assimilation for distributed hydrological catchment modeling via ensemble Kalman filter. Adv. Water Resour. 33, 678–690. Xu, C.Y., 1999. Estimation of parameters of a conceptual water balance model for ungauged catchments. Water Resour. Manag. 13, 353–368. Young, P.C., 2002. Advances in real-time flood forecasting. Phil. Trans. R. Soc. Lond. A 360, 1433–1450. Yu, Z., 2000. Assessing the response of subgrid hydrologic processes to atmospheric forcing with a hydrologic model system. Global Planet Change 25, 1–17. Yu, Z., Lü, H., Zhu, Y., Drake, S., Liang, C., 2010. Long-term effects of revegetation on soil hydrological processes in vegetation-stabilized desert ecosystems. Hydrol. Process. 24, 87–95. Zhao, R.J., 1992. The Xinanjiang model applied in China. J. Hydrol. 135, 371–381. Zhao, R.J., Zhang, Y.L., Fang, L.R., Liu, X.R., Zhang, Q.S., 1980. The Xinanjiang model. In: Hydrological Forecasting Proceedings Oxford Symposium. IAHS. pp. 351– 356.