Uncertainty and its propagation estimation for an integrated water system model: An experiment from water quantity to quality simulations

Journal of Hydrology 565 (2018) 623–635

Contents lists available at ScienceDirect

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

Research papers

Uncertainty and its propagation estimation for an integrated water system model: An experiment from water quantity to quality simulations

T

⁎

Yongyong Zhanga, , Quanxi Shaob a

Key Laboratory of Water Cycle and Related Land Surface Processes, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China b CSIRO, Data 61, Private Bag 5, Wembley, Western Australia 6913, Australia

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by G. Syme, Editor-in-Chief, with the assistance of Jesús Mateo-Lázaro, Associate Editor

Multiple uncertainty sources directly cause inaccurate simulations for water related processes in complicated integrated models, as such models include many interactive modules. A majority of existing studies focus on the uncertainties of parameter and model structure, and their effects on the model performance for a single process (e.g., hydrological cycle or water quality). However, comprehensive uncertainties of different modules and their propagations are poorly understood, particularly for the integrated water system model. This study proposes a framework of uncertainty and its propagation estimation for integrated water system model (HEQM) by coupling the Bootstrap resampling method and SCE-UA auto-calibration technique. Parameter and structure uncertainties of both hydrological cycle and water quality modules are estimated, including final distributions of parameters and simulation uncertainty intervals. Additionally, the effect of uncertainty propagation of hydrological parameters is investigated. Results show that: (1) HEQM simulates daily hydrograph very well with the coefficient of efficiency of 0.81, and also simulates the daily concentrations of ammonia nitrogen satisfactorily with the coefficient of efficiency of 0.50 by auto-calibration in the case study area; (2) The final ranges of all interested hydrological parameters are reduced obviously, and all the parameter distributions are well-defined and show skew. The uncertainty intervals of runoff simulation at the 95% confidence level bracket 18.7% of all the runoff observations due to uncertainties of parameter, and 86.0% due to both parameter and module structure, respectively; (3) The uncertainty propagation of hydrological parameters changes the optimal values of 37.5% of interested water quality parameters, but does not obviously change the water quality simulations which match well with the prior simulations throughout the period and bracket only 1.7% of observations at the 95% confidence level. Due to the further introduction of module structure uncertainties, 94.8% of observations are bracketed, only except the extreme high and low water quality concentrations; (4) The uncertainty of water quality parameters contributes 12.1% of total water quality simulations at the 95% confidence level. The figure increases to 21.0% and 92.0% if the uncertainty propagation of hydrological parameters, structure uncertainties of water quality module are considered, respectively. Therefore, although the parameter uncertainty and its propagation contribute a certain proportion of the whole simulation uncertainties, the module structure itself is the primary uncertainty source for the integrated water system model (HEQM), particular for the water quality modules.

Keywords: Uncertainty propagation Model performance Bootstrap resampling SCE-UA auto-calibration HEQM

1. Introduction Linkages, interconnections and interdependencies of water cycle have been gradually recognized at basin or global scales due to the rapid growth of environmental science and the constantly emerging water issues (e.g., drought, flooding, erosion, pollution and ecological degradation) (GWSP, 2005). Integrated consideration or simulation of multiple water related processes become a new trend along with further

⁎

explorations of interaction mechanisms among multiple processes, rapid developments of computer facilities and observation techniques of multiple data sources (Paola et al., 2006; Zhang et al., 2016a,b). Many successful model integrations have been implemented with different objectives in the earth system studies. For example, land surface models (e.g., VIC-Variable Infiltration Capacity; Liang et al., 1994) could be coupled with hydrological models to reveal interactions and feedbacks between atmosphere and hydrology at large scale. Similarly,

Corresponding author. E-mail address: [email protected] (Y. Zhang).

https://doi.org/10.1016/j.jhydrol.2018.08.070 Received 22 May 2018; Received in revised form 27 August 2018; Accepted 28 August 2018 Available online 31 August 2018 0022-1694/ Crown Copyright © 2018 Published by Elsevier B.V. All rights reserved.

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

1992; Bates and Campbell, 2001; Engeland and Gottschalk, 2002; Montanari and Brath, 2004; Li et al., 2010a,b; Shao et al., 2014; Arsenault et al., 2015) and water quality models (Beck, 1987; Freni et al., 2008; Jia et al., 2018; Novic et al., 2018). The comprehensive uncertainty estimations in the integrated model of multiple processes are still deficient. Furthermore, most of existing studies about uncertainty propagation analysis are limited to the impact assessments of observation quality on model performances which are induced by potential errors from sampling, instrument, laboratory, observation method or algorithm (Harmel et al., 2006; Xu et al., 2006; Leta et al., 2015). The related uncertainties include uncertainties associated with model inputs (e.g., climate variables, geographic information) (Crosetto et al., 2001; Gabellani et al., 2007; Shao et al., 2012; Yen et al., 2014; Novic et al., 2018) and observations used for the model calibration (e.g., flow regime, water quality variables) (Harmel et al., 2006; Shao et al., 2014; Yen et al., 2014). Both of two uncertainty sources are propagated and probably distort the probability distributions of model parameters, and thus disturb the capability of integrated models to portray the real world, particularly for the simulations of subsequent modules. For example, Harmel et al. (2006) examined that the probable uncertainties for observed variables ranged from 6% to 19% for streamflow, from 11% to 100% for NH4-N, from 11% to 104% for total nitrogen (TN), and from 8% to 110% for total phosphorous (TP), among which the contributions of sample collection, preservation storage and laboratory analysis were from 4% to 48%, from 2% to 16% and from 5% to 21%, respectively. However, only a few studies are reported about the uncertainty propagation investigation among different modules, particularly for integrated water system models. The main purpose of this study is to comprehensively assess the uncertainties of multiple modules and their propagations among different modules for complicated integrated water system models. As a typical integrated water system model, HEQM (Hydrological, Ecological and water Quality Model) is adopted to investigate the effect of multiple uncertainty sources on parameter distributions and simulation performances of both hydrological and water quality modules. The specific objectives are to: (1) propose a comprehensive assessment of uncertainty sources and their propagations for HEQM by using Bootstrap resampling with SCE-UA optimization technique; (2) estimate the probability distributions of hydrological parameters and uncertainty intervals of runoff simulation caused by uncertainties of parameters and module structures; (3) estimate the probability distributions of water quality parameters and uncertainty intervals of water quality simulation caused by parameter uncertainty propagation of hydrological cycle module; (4) estimate the probability distributions of water quality parameters and uncertainty intervals of water quality simulation caused by uncertainties of parameter and module structure. This study is expected to extend the scope of model uncertainty analysis, and assist modelers in further improvements and calibrations of complicated integrated models of multiple processes.

hydrological models (e.g., SWAT-Soil and Water Assessment Tool; Arnold et al., 1998) could be coupled with soil erosion or biogeochemical processes to reveal precipitation-induced losses of water, soil and nutrient from lands to river networks. Hydrological models are also able to be couple with hydrodynamic and water quality models of water bodies (i.e., rivers or lakes) to capture the migrations of waterborne variables with high spatial and temporal resolutions, such as EFDC (Environmental Fluid Dynamic Code) (Hamrick, 1992). Ecosystem models (EPIC- Erosion/Productivity Impact Calculator; Sharpley and Williams, 1990; DNDC: DeNitrification/DeComposition; Li et al., 1992) could also be coupled with evapotranspiration model and soil biogeochemical model to reveal vegetation growth processes with considerations of nutrient and water stresses. Furthermore, integrated water system model is proposed and usually formed as a chain of water related modules to capture the interactions and feedbacks among physical, biological and geochemical processes, as well as the impacts of water-related human activities (GWSP, 2005), such as CLM series (Community Land Model) (Dai et al., 2003). However, not all the processes are physically interpreted by mathematical equations due to the current insufficient knowledge (Willems, 2008). Empirical conceptualization and theoretical simplification are usually adopted and are easy to introduce potential uncertainties of parameters and model structures into the complicated models (Todini, 2007; Freni et al., 2009). For the integrated models that are made up of many modules, multiple uncertainty sources of upstream modules are transferred to the downstream modules as inputs (Freni et al., 2008). Therefore, along with the increasing of coupled modules, multiple uncertainty sources from different modules not only affect the simulation performance of their own modules, but also might be accumulated and propagated to subsequent modules and thus distort their performances and even the whole model. For example, in the applications of integrated water quantity and quality models, water quality simulation performance is not usually satisfying using the step-by-step calibration approach even though the water quality modules are well formulated, which is probably caused by the uncertainty or error propagation from upstream modules (Zhang et al., 2016b). It is critical to investigate the estimations of uncertainty sources of different modules and their propagation, as well as their effects on simulation performance. Large majorities of studies are implemented to identify model uncertainty sources, and assess their effects on model performance. The identified uncertainty sources mainly include model parameter uncertainty (Beven and Binley, 1992; Bates and Campbell, 2001; Beven and Freer, 2001; Yang et al., 2007; Li et al., 2010a,b; Leta et al., 2015) and model structure uncertainty (Beven and Binley, 1992; Refsgaard et al., 2006; Li et al., 2010a). Most of the existing techniques are categorized into two classes, i.e., (1) the frequentist approach with model calibration techniques (e.g., Shuffled Complex Evolution: SCE-UA, Particle Swarm Optimization: PSO) which is advantageous to be implemented without timing consumption and the representative techniques are GLUE procedure (Generalized Likelihood Uncertainty Estimation) (Beven and Binley, 1992; Beven and Freer, 2001) and Parameter Solution (Duan et al., 1992), both of which are based on subjective determination of generalized likelihood measures between simulations and observations, Sequential Uncertainty Fltting algorithm (SUFI-2) with global sampling techniques based on multi-criteria thresholds of model calibration (Abbaspour et al., 2007), and Bootstrap resampling with recalibration (Li et al., 2010b; Novic et al., 2018); (2) the classical Bayesian theorem with sampling techniques (e.g., Markov Chain Monte Carlo, Latin hypercube) based on the observations and prior information of model parameters (Bates and Campbell, 2001; Engeland and Gottschalk, 2002; Montanari and Brath, 2004; Yang et al., 2007), which is robust and widely used to estimate the reliable uncertainties of model parameters. However, all of these studies are model-specific (Engeland and Gottschalk, 2002; Gallagher and Doherty, 2007), and only focus on the inherent uncertainties of single process model or model structure, e.g., hydrological models (Beven and Binley,

2. Models and methodology 2.1. Integrated water system model (HEQM) HEQM is an integrated water system model proposed by Zhang et al. (2016a) in order to investigate hydrological cycle processes, its accompanied biogeochemical and water quality processes as well as their interactions at catchment scale. The main water related processes are mathematically described by hydrological cycle module (HCM), soil erosion module (SEM), overland water quality module (OQM), water quality module in water bodies (WQM), crop growth module (CGM), soil biochemical module (SBM) and dams regulation module (DRM) (Fig. S1 in the Supplementary material). Furthermore, a parameter analysis tool (PAT) is provided to conveniently conduct the parameter sensitivity analysis, model calibration and performance assessment. All 624

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

−0.10 and 0.10 for Nbias, greater than 0.85 for r and greater than 0.75 for NS. For the NH4-N concentration simulation, the performance is at the satisfactory rating if the statistic values are between 0.30 and 0.50 or between −0.50 and −0.30 for Nbias, between 0.75 and 0.80 for r and between 0.50 and 0.60 for NS, and the performance is at the very good rating if the statistic values are between −0.15 and 0.15 for Nbias, greater than 0.85 for r and greater than 0.75 for NS.

the equations are given in Zhang et al. (2016a), and the detailed interactions between hydrological cycle and water quality processes are summarized in Zhang et al. (2016b) and the Supplementary material. In order to consider spatial heterogeneities of catchment attributes and vegetation coverage, three categories of spatial calculation units are designed, including sub-basin/catchment, land use/land cover and crop. The sub-basin/catchment unit is divided based on the digital elevation model (DEM) using the spatial analyst tools of ArcGIS software. The sub-basin/catchment unit is divided into the land use/land cover unit according to the national land use/land cover classification (China’s national standard, 2007) and seven user-defined types of land use/land cover are considered for each sub-basin/catchment. The land use/land cover unit (e.g., paddy, dryland agriculture, orchard) where the agricultural activities occur, is further divided into the crop unit. The detailed functions of different spatial calculation units and their related modules are given in Zhang et al. (2016a).

2.2.2. Auto-optimization algorithm and Bootstrap method Shuffled Complex Evolution method developed at the University of Arizona (SCE-UA) is a powerful global optimization strategy based on a synthesis of best features of controlled random search, competitive evolution and coupled with a newly complex shuffling (Duan et al., 1994). This strategy is originally designed to solve the auto-optimization problems of water related models (Duan et al., 1992), and has been widely accepted so far due to its efficiency and effectiveness, as well as its open sources (Arsenault et al., 2015; Fowler et al., 2016). Bootstrap is a nonparametric method for statistic estimation and its distribution testing proposed by Professor Efron (1979). The major feature of this method is that the sampling distribution of some prespecified random variables can be estimated on the basis of the original samples, even though the true probability distribution is unknown. It is a sample augmentation method and able to fully utilize the information of original samples, which is strongly advantageous to solve the sampling problems when only small samples are available. It has been widely used in hydrological community, such as uncertainty analysis (Li et al., 2010b; Shao et al., 2012), time series analysis (e.g., trend, correlation, frequency) (Yue and Wang, 2002; Hirsch et al., 2015), hydrological event design (e.g., storm or flood) (Zucchini and Adamson, 1989; Hu et al., 2015) and risk analysis (Wucchini and Adamson, 1988). For both SCE-UA and Bootstrap methods, sufficient resampling is needed and the specific resampling of both two methods in our study is introduced in Section 2.3.

2.2. Framework of uncertainty propagation analysis 2.2.1. Model residual and model performance assessment Model residual is the discrepancies between the simulated and observed variables. It is formulated as

⎧ ε (t ) = O (t )−S (t ) , ⎨ S (t ) = f (θ )̂ ⎩

(1)

whereε (t ) is the simulated variable residual; O(t) and S(t) are the observed and simulated variables at time t. The interested simulated variables are runoff and ammonia nitrogen concentration (NH4-N) in our study. f (θ )̂ is the mathematical expression for the variables from relevant HEQM modules with θ ̂ being the estimated value of unknown q q q q parameter θ . For runoff simulation, θ ̂ = (θ ̂ 1, θ ̂ 2, ⋯, θ ̂ pq) and pq is the number of interested hydrological parameters. For the NH4-N sic c c c mulation, θ ̂ = (θ ̂ 1, θ ̂ 2, ⋯, θ ̂ pc ) and pc is the number of interested water quality parameters. Root mean square error (RMSE) is adopted as the objective function for the auto-optimizations of hydrological and water quality parameters because it is straightforward and widely used (Zhang et al., 2016b). The model performance is close to the optimal if RMSE approaches to 0.0. Additionally, three other statistics are used to further evaluate the model performances including normalized bias (Nbias), correlation coefficient (r) and coefficient of efficiency (NS) because several simulation performance ratings are usually divided based on these statistics, such as very good, good, satisfactory and unsatisfactory according to the existing studies (e.g., Moriasi et al., 2007; Zhang et al., 2016b). All the equations are given as follows.

Root mean square error(RMSE ): RMSE =

N

∑i =1 (Oi−Si)2

N

2.2.3. Integrated framework In our study, three major uncertainty sources from water quantity and quality simulations are considered for the uncertainty analysis estimation of HEQM, i.e., uncertainties for hydrological cycle modules, parameter uncertainty propagation analysis and uncertainties of water quality modules (Fig. 1). The SCE-UA algorithm has been coupled in the PAT of HEQM, which is able to achieve the value exchanges of parameters and objectives online. Additionally, the coupling between Bootstrap method and HEQM is off-line by reading in and out the resampled residual groups in the text file. The Bootstrap resampling is conducted using the sample function of base Package (version 3.3.0) (Becker et al., 1988) in the R platform (version 3.1.1) (R Development Core Team, 2010).

(2) 2.2.3.1. Uncertainty analysis for hydrological cycle modules. In this part, hydrological parameter estimations are conducted, and uncertainties of parameters and module structures are analyzed for the effects on runoff simulation performance. The detailed procedure is given as follows.

N

Normalized bias(Nbias ): Nbias =

Correlation coefficient(r ): r =

∑i = 1 (Oi−Si ) N

∑i = 1 Oi

(3)

N ∑i = 1 (Oi−O¯ )·(Si−S¯ ) N N ∑i = 1 (Oi−O¯ )2 · ∑i = 1 (Si−S¯ )2

(1) Select the interested parameter groups θ q = (θ q1, θ q2, ⋯, θ q pq) of hydrological cycle module (HCM) (Table 1); q (2) Obtain the optimal values of hydrological parameter groups θ ̂ by the SCE-UA auto-calibration, and their corresponding optimal series of runoff Sq(t) (t = 1,2,…,N); (3) Calculate the simulation residuals of runoff ε q (t ) using Eq. (1), and resample the residuals with replacement for maxb times by the Bootstrap method to form new residual series εnq (t , i) (i = 1,2,…, maxb); (4) Add εnq (t , i) to the optimal runoff series and form the new series Snq (t , i) = εnq (t , i) + S q (t ) (i = 1,2,…, maxb); (5) Calibrate the HEQM by SCE-UA auto-optimization again according

(4)

N

Coefficient of efficiency(NS ): NS = 1−

∑i = 1 (Oi−Si )2 N ∑i = 1 (Oi−O¯ )2

(5)

where O¯ and S¯ are the average values of observed and simulated variables, respectively; N is the total length of time series. According to the model performance evaluation in Zhang et al. (2016b), the runoff simulation performance is at the satisfactory rating if the statistic values are between 0.15 and 0.25 or between −0.25 and −0.15 for Nbias, between 0.75 and 0.80 for r, and between 0.50 and 0.60 for NS, and the performance is at the very good rating if the statistic values are between 625

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Fig. 1. Schematic diagram of uncertainty analysis and its propagation by coupling HEQM with Bootstrap resampling and SCE-UA auto-calibration.

Sc(t) (t = 1,2,…,N), and the simulation residualsε c (t ) ; (3) Re-calibrate the interested water quality groups θ c = (θ c1, θ c 2, ⋯, θ c pc ) by fixing the hydrological parameters q θ q = (θ q1, θ q2, ⋯, θ q pq) as the optimal values θn̂ (i) in Step (5) of

to the new series Snq (t , i) (i = 1,2,…,maxb) and obtain the correq sponding optimal values of hydrological parameter groups θn̂ (i) q and simulated runoff series S ̂ n (t , i) (i = 1,2,…, maxb); (6) Obtain the marginal distributions of interested hydrological paraq meters θn̂ and the uncertainty intervals of runoff by the two-side q q q confidence intervals at level α, i.e., Spq = [Sn̂ , α /2 (t ), Sn̂ ,1 − α /2 (t )]; (7) Resample the residuals with replacement for maxb times to form maxb groups of new residual series εmq (t , i) again, and repeat Step q q (4) to form new runoff values S ̂ m (t , i) = S ̂ n (t , i) + εmq (t , i) (i = 1,2,…, maxb); (8) Obtain the uncertainty intervals of runoff further by module structure and parameter uncertainties at level α, i.e., q q q ̂ ̂ Spq + mq = [Sm, α /2 (t ), Sm,1 − α /2 (t )];

Section 2.2.3.1 individually for maxb times, and obtain the optimal c values of water quality parameter groups θn̂ (i) (i = 1,2,…, maxb) c and their corresponding water quality series S ̂ n (t , i) (i = 1,2, …,maxb); c (4) Obtain the marginal distributions of water quality parameters θn̂ , and the uncertainty propagation intervals of water quality variables due to hydrological parameter uncertainties by the two-side conc c c _opt ̂ ̂ fidence interval at level α, i.e., Spc , pq = [Sn, α /2 (t ), Sn,1 − α /2 (t )]; (5) Resample the residuals ε c (t ) with replacement for maxb times to form new residual series εnc (t , i) (i = 1,2,…, maxb) for water quality variables, and form new values of water quality variables c c S ̂ m (t , i) = S ̂ n (t , i) + εnc (t , i) (i = 1,2,…, maxb); (6) Obtain the uncertainty intervals of water quality variables due to parameter uncertainty propagation of hydrological cycle modules, parameter and structure uncertainty of water quality modules at c c c _opt ̂ ̂ level α, i.e., Spc + mc, pq = [Sm, α /2 (t ), Sm,1 − α /2 (t )].

2.2.3.2. Parameter uncertainty propagation analysis for hydrological cycle modules. In this part, parameter uncertainty propagation of hydrological cycle modules is analyzed for the effects on distributions of water quality parameters, and water quality simulation performance. The detailed procedure is given as follows. (1) Select the interested parameter groups θ c = (θ c1, θ c 2, ⋯, θ c pc ) of water quality module (WQM) (Table 1); (2) Auto-calibrate the water quality parameters θ c = (θ c1, θ c 2, ⋯, θ c pc ) by fixing the hydrological parameters θ q = (θ q1, θ q2, ⋯, θ q pq) as the q optimal values θ ̂ , and obtain the optimal values of water quality c parameter groups θ ̂ , their corresponding optimal series of runoff

2.2.3.3. Uncertainty analysis for water quality modules. In this part, water quality parameter estimations are conducted, including parameter uncertainty, module structure uncertainty, and their effects on performance of water quality simulation. The detailed procedure is given as follows.

Table 1 Selected sensitive hydrological and water quality parameters of HEQM for autocalibration and uncertainty analysis. Modules

Name

Units

Min

Max

Definition

Hydrological cycle modules

Wfc Wsat g1 g2 KET Kr Tg Kg fc Kset(NH4) Rd(COD) Kd(NH4) Rset(COD) Kd(COD) Rd(NH4) Rset(NH4) Rset(orgN)

none none none none none none days none mm/hour m/year 1/day 1/day 1/day 1/day 1/day m/year 1/day

0.200 0.450 0.000 0.000 0.000 0.000 1.00 0.000 0.00 0.00 0.020 0.100 −0.360 0.020 0.100 −50.00 0.001

0.450 0.750 3.000 3.000 3.000 1.000 100.00 1.000 120.00 100.00 3.400 1.00 0.360 3.400 1.000 50.00 0.1

Field capacity of soil Saturation moisture capacity of soil Basic runoff coefficients for different land use types Influence coefficient of soil moisture for different land use types Adjustment factor of evapotranspiration Interflow yield coefficient Delay time for aquifer recharge Baseflow yield coefficient Steady state infiltration rate of soil Settling rate of NH4-N at 20 °C in reservoir COD deoxygenation rate at 20 °C in channel Bio-oxidation rate of NH4-N at 20 °C in reservoir COD settling rate at 20 °C in channel COD deoxygenation rate at 20 °C in reservoir Bio-oxidation rate of NH4-N at 20 °C in channel NH4-N settling rate at 20 °C in channel Organic N settling rate at 20 °C in channel

Water quality modules

626

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

(1) Auto-calibrate the water quality parameters θ c = (θ c1, θ c 2, ⋯, θ c pc ) by fixing the hydrological parameters θ q = (θ q1, θ q2, ⋯, θ q pq) as the q optimal values θ ̂ , and obtain the optimal values of water quality c parameter groups θ ̂ , and their corresponding optimal series of

mainstream (Zhang et al., 2015). The detailed procedure of HEQM setup for Fuyang controlled catchment is introduced in Zhang et al. (2016b), involving sub-catchment delineation, and interpolation of climate and socio-economic data. The essential data sources are also provided in the Supplementary material. Paddy, dryland agriculture, forest, grass, water, urban and unused lands are considered as the land use/land cover units of HEQM. In our study, the daily runoff and NH4-N concentration from 2003 to 2005 are selected as the simulated variables of water quantity and quality, respectively, and the first half year (from 1st January 2003 to 1st June 2003) is set as the warm up period for water quality simulation. In the integrated framework, nine interested hydrological parameters and eight water quality parameters of HEQM (Table 1) are selected according to the previous studies (e.g., Zhang et al., 2016b) for calibration and uncertainty analysis, i.e., pq = 9 and pc = 8. The related parameters are soil field capacity (Wfc), soil saturation moisture capacity (Wsat), basic runoff coefficient (g1) and influence coefficient of soil moisture (g2) for different land use types, adjustment factor of evapotranspiration (KET), interflow yield coefficient (Kr), delay time for aquifer recharge (Tg), baseflow yield coefficient (Kg), soil steady state infiltration rate (fc) for the hydrological cycle modules, and settling rate of NH4-N at 20 °C in reservoir (Kset(NH4)), COD deoxygenation rate at 20 °C in channel (Rd(COD)), bio-oxidation rate of NH4-N at 20 °C in reservoir (Kd(NH4)), COD settling rate at 20 °C in channel(Rset(COD)), COD deoxygenation rate at 20 °C in reservoir (Kd(COD)), bio-oxidation rate of NH4-N at 20 °C in channel(Rd(NH4)), NH4-N settling rate at 20 °C in channel(Rset(NH4)), organic N settling rate at 20 °C in channel (Rset(orgN)) for the water quality modules. Furthermore, 1000 samplings of both runoff and water quality simulation residual groups are implemented by the Bootstrap method and the significant level is set as 5%, i.e., maxb = 1000 and α = 5%. In the SCE-UA algorithm, the maximum iterations are set as 2000 for runoff calibration and 6000 for water quality calibration, respectively, while the other parameters are set as the default values. In order to enhance the calculation efficiency, 20 parallel optimizations are implemented with 50 residual groups each.

runoff Sc(t) (t = 1,2,…,N); (2) Calculate the simulation residuals of water qualityε c (t ) using Eq. (1), and resample the residuals with replacement for maxb times by Bootstrap method to form new residual series εnc (t , i) (i = 1,2,…, maxb); (3) Add εnc (t , i) to the optimal water quality series and form the new series Snc (t , i) = εnc (t , i) + S c (t ) (i = 1,2,…, maxb); (4) Calibrate the HEQM by SCE-UA auto-optimization according to the new series Snc (t , i) (i = 1,2,…, maxb) and obtain the corresponding c optimal values of water quality parameter groups θn̂ (i) and water c quality variable series S ̂ n (t , i) (i = 1,2,…, maxb); (5) Obtain the marginal distributions of interested water quality parac meters θn̂ , and the uncertainty intervals of water quality variable simulation by the two-side confidence intervals at level α, i.e., c c c Spc = [Sn̂ , α /2 (t ), Sn̂ ,1 − α /2 (t )]; (6) Resample the residuals with replacement for maxb times to form new residual series εmc (t , i) again, and repeat Step (4) to form new c c water quality values S ̂ m (t , i) = S ̂ n (t , i) + εmc (t , i) (i = 1,2,…, maxb); (7) Obtain the uncertainty intervals of water quality further by module structure and parameter uncertainties at level α, i.e., c c c ̂ ̂ Spc + mc = [Sm, α /2 (t ), Sm,1 − α /2 (t )]. 2.2.4. Assessment of model uncertainty intervals Besides the statistics (Nbias, r and NS) of model performance assessment mentioned-above, the other two statistics presented by Li et al. (2010a) are adopted to assess the uncertainty intervals of both runoff and water quality simulations, e.g., the average relative interval length (ARIL) at the 95% confidence level, and the percent of observations bracketed by the 95% confidence interval (P-95CI). The equations are given as follows.

ARIL =

1 N

P−95CI =

N

∑t=1 (S97.5% (t )−S2.5% (t ))/O (t )

(6)

NOin × 100% N

(7)

3. Results 3.1. Model auto-calibration

where S97.5% and S2.5% are the 97.5th and the 2.5th percentages of maxb groups of runoff or water quality simulation series; NOin is the number of observations which are bracketed in the 95% confidence interval. Smaller values of ARIL indicate the uncertainty intervals are very narrow, and the greater values of P-95CI indicate the uncertainty intervals are very reliable (Jin et al., 2010; Li et al., 2010b). The optimal values of ARIL and P-95CI are 0.00 and 95%, respectively.

The selected hydrological and water quality parameters of HEQM are auto-calibrated step-by-step using the SCE-UA approach. The selected hydrological (θ q ) and water quality (θ c ) parameters are conq c verged to the optimal values, i.e., θ ̂ and θ ̂ . The simulated daily hydrograph matches very well with the observations with the Nbias, r and NS of 0.01, 0.90 and 0.81 (Fig. 2a and Table 2). The simulated NH4-N performance is relatively good for the whole records with the Nbias, r and NS of 0.26, 0.80 and 0.50, respectively. The major deviations between simulations and observations are in the pre-flood season (January to May), when the simulated NH4-N concentrations are underestimated (Fig. 2b). The probable explanation is that HEQM neglects river sediment resuspension, which is a frequent event and contributes high pollutant concentration during the pre-flood season due to the sluice regulation for the flood discharge (Zhang et al., 2015). According to the evaluation criteria in Moriasi et al. (2007) and Zhang et al. (2016b), the HEQM performance is at the very good rating for daily runoff simulation and at the satisfactory rating for daily NH4-N concentration simulation, respectively.

2.3. Study area and model setup Shaying River catchment (111°56′–116°31′E, 32°28′–34°54′N) is the largest catchment of Huai River basin with the total drainage area of 36,651 km2 (Fig. S2 in the Supplementary material). As the major grain producing area of the Huai River basin, the dominating land use type is the agricultural land accounting for 85% of the total catchment area. It is also the congeries of highly polluted cottage industries such as papermaking, fertilizer, leather, chemistry and dipdye. Over 0.25 tons of industry wastewater per year are discharged into Shaying River, in which the emission loads of NH4-N and chemical oxygen demand (COD) are beyond the water environmental capacity of most streams (Zhang et al., 2015). Fuyang station (115°50′E, 32°54′N), located at the downstream of Shaying River, is selected as the outlet of our study area. Its control area accounts for 94.9% of the whole Shaying River catchment. It is a critical hydrological and water quality station of Shaying River, where water quality condition obviously affects that of the Huai

3.2. Uncertainty estimation of hydrological cycle modules 3.2.1. Hydrological parameter estimations The marginal final distributions of all the selected hydrological 627

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

The runoff simulation uncertainties due to the uncertainties of parameter and module structure are also considered (Fig. 4b). The simulation uncertainty intervals at the 95% confidence level (i.e., ARIL = 3.66) are much wider than those intervals only due to the parameter uncertainties. Most of observations are bracketed (i.e., P95CI = 86.0%) (Table 4), except the peak flow observations in the flood season due to their underestimations by HEQM. Therefore, the improvement of hydrological cycle module is likely to capture more runoff observations, but is still disadvantageous to predict the high flow whose uncertainty is probably caused by the uncertainties of model inputs or runoff observations, as well as the regulation of downstream reservoirs and sluices. Nonetheless, the statistics of model performance assessment are from −0.10 to 0.04 for Nbias, from 0.84 to 0.89 for r, and from 0.70 to 0.78 for NS. The runoff simulation performance of uncertainty intervals is still good. 3.3. Propagation effects of hydrological parameter uncertainties on water quality simulation 3.3.1. Effects on estimations of water quality parameters All the selected water quality parameter groups are auto-calibrated by fixing the hydrological parameter groups at the optimal values which are calibrated according to the new runoff series. The marginal final distributions of water quality parameters are given in Fig. 5. Compared to the initial ranges of all the selected parameters, the ranges of final distributions are reduced remarkably with the ratios from 30.7% to 100.0%, and the reductions of Kset(NH4), Rd(COD), Kd(NH4) and Rd(NH4) are the most remarkable. The final distributions of selected parameters also obey the skew distributions. Except Kd(COD), the distributions of all the water quality parameters show the positive skewness with the values from 0.774 (Rd(COD)) to 9.254 (Kset(NH4)). For Kset(NH4), Kd(NH4), Rset(COD), Rd(NH4) and Rset(NH4), most of their optimal values according to the new runoff series are very close to their initial optimal values. However, most of the optimal values of Rd (COD), Rset(orgN) are close to their initial minimum values and only that of Kd(COD) is close to the initial maximum value. Both of the minimum and maximum values are quite different from their initial optimal values. Therefore, the uncertainty propagation of hydrological parameters does not seem to affect the distributions of Kset(NH4), Kd (NH4), Rset(COD), Rd(NH4) and Rset(NH4), but obviously affect the distributions of Rd(COD), Rset(orgN) and Kd(COD).

Fig. 2. Comparison of simulations and observations for the daily runoff (a) and water quality concentration (b) from 2013 to 2017 at Fuyang station. Table 2 Performance assessment of runoff and water quality simulations. Criteria Nbias r NS

Runoff *

0.01 0.90* 0.81*

Water quality 0.26 0.80** 0.50***

Note: The values with superscripts “*”, “**”, and “***” are the simulation performance at the very good, good and satisfactory ratings, respectively.

3.3.2. Effects on water quality simulation performance The uncertainty intervals of water quality simulation at the 95% confidence level are quite narrow with the ARIL of 0.045 (Fig. 6a), and are very close to the optimal simulated values throughout the period, except in April and May. Only 1.7% of observations are bracketed in the uncertainty intervals (i.e., P-95CI = 1.7%). Moreover, the criteria of simulation performance assessment are from 0.17 to 0.23 for Nbias, from 0.75 to 0.78 for r, and from 0.48 to 0.55 for NS (Table 4). Compared with the initial auto-calibration results, the Nbias is improved slightly while the r is weakened slightly. Therefore, the uncertainty propagation of hydrological parameters also does not affect the performance of water quality simulation obviously. Furthermore, Fig. 6b presents the uncertainty intervals of water quality simulations due to the uncertainty propagation of hydrological cycle parameters and uncertainty of water quality module structure. With the introduction of module structure uncertainty, the ARIL increases a lot (i.e., 5.94), which is much wider than the intervals only caused by the parameter uncertainty propagation. A large majority of observations (P-95CI = 94.8%) are bracketed in the uncertainty intervals, which is very close to the optimum of P-95CI (95%). Thus, the uncertainty intervals of water quality simulations are very robust by our proposed integrated framework. However, the extremely high and low observations of water quality concentrations are still not captured and the probable explanations are the uncertainties of model inputs or

parameters are well-defined and obey the skew distributions (Fig. 3), among which those of WMc, WM, Kg2, Kr, Tg and Krg show the positive skewness with the values from 0.061 (WMc) to 1.023 (Kg2), and the others show the negative skewness with the values from −0.022 (fc) to 0.530 (KETp) (Table 3). Compared with the initial ranges of parameter values, all the final ranges at the 95% confidence interval are shortened remarkably. The change of KETp is the greatest with the reduction ratio of 70.1% while the change of Krg is the least with the reduction ratio of 19.0%. Additionally, the optimal parameter value is close to the median value for WMc, Kg1, KETp, Kr, Tg and fc, to the minimum value for WM, and to the maximum value for Kg2 and Krg.

3.2.2. Uncertainty intervals of runoff simulation Fig. 4a shows the comparison between runoff observations and simulations at the 95% confidence interval due to parameter uncertainties. The ARIL is 0.37, which is close to the corresponding optimal value (0.0). It indicates that the uncertainty intervals of runoff simulation due to parameter uncertainties are quite narrow. Furthermore, although only a few observations are bracketed in the uncertainty intervals (i.e., P-95CI = 18.7%), the simulation performance is still very good with the Nbias, r and NS from −0.08 to 0.01, 0.89 to 0.90, and 0.79 to 0.81, respectively. 628

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Fig. 3. Marginal posterior densities of all the selected hydrological parameters of HEQM.

3.4. Uncertainty estimation of water quality modules

runoff observations, as well as the module structures. Therefore, the existing water quality modules of HEQM still have a large potential to capture most of the observations by improving their module structures, e.g., adopting more physical formulas to replace the conceptual or simplified formulas.

3.4.1. Water quality parameter estimations The marginal final distributions of all the selected water quality parameters are well-defined and obey the skew distributions (Fig. 7). The final distributions of Kset(NH4), Kd(NH4), Kd(COD), Rd(NH4) and Rset(orgN) show the positive skewness with the values from 0.102 to 11.930, and the others (i.e., Rd(COD), Rset(COD), Rset(NH4)) show the

629

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Table 3 Optimal values and posterior parameter distributions of all the selected hydrological and water quality parameters of HEQM. Parameters

Hydrological parameters

Water quality parameters due to uncertainty propagation of hydrological parameters

Water quality parameters

Name

Wfc Wsat g1 g2 KET Kr Tg Kg fc Kset(NH4) Rd(COD) Kd(NH4) Rset(COD) Kd(COD) Rd(NH4) Rset(NH4) Rset(orgN) Kset(NH4) Rd(COD) Kd(NH4) Rset(COD) Kd(COD) Rd(NH4) Rset(NH4) Rset(orgN)

Optimum

0.307 0.450 2.002 2.991 0.444 0.265 66.34 0.999 45.93 0.000 0.022 0.100 −0.307 1.671 0.100 −50.00 0.030 0.000 0.022 0.100 −0.307 1.671 0.100 −50.00 0.030

Posterior parameter distribution Min

2.5th Percentile

97.5th Percentile

Max

Variance

Skewness

Reduction of range

0.201 0.450 0.359 0.071 0.004 0.180 1.64 0.217 7.01 0.000 0.020 0.100 −0.360 0.896 0.100 −50.00 0.001 0.000 0.024 0.100 −0.360 0.021 0.100 −49.99 0.001

0.223 0.466 0.755 0.482 0.044 0.227 21.07 0.306 19.15 0.000 0.021 0.100 −0.360 1.056 0.100 −50.00 0.001 0.001 0.149 0.100 −0.349 0.123 0.100 −49.79 0.001

0.405 0.697 2.865 2.365 0.725 0.429 83.69 0.939 112.21 0.020 0.044 0.102 −0.307 3.398 0.102 −45.67 0.040 0.745 3.397 0.102 0.315 3.214 0.101 49.92 0.051

0.445 0.748 2.999 2.991 0.830 0.823 99.34 0.999 119.84 0.157 0.076 0.109 −0.307 3.400 0.110 −37.59 0.097 1.170 3.400 0.120 0.358 3.390 0.107 50.00 0.077

0.046 0.058 0.530 0.451 0.179 0.054 15.912 0.170 24.277 0.008 0.007 0.001 0.010 0.872 0.001 1.318 0.012 0.215 1.164 0.001 0.175 0.843 0.000 39.428 0.012

0.061 0.360 −0.274 1.023 −0.530 0.907 0.202 0.361 −0.022 9.254 0.774 4.957 4.805 −1.217 6.362 3.942 2.572 1.301 −1.287 12.816 −0.175 0.102 11.930 −0.765 2.114

27.3% 22.9% 29.7% 37.2% 77.3% 79.8% 36.8% 36.7% 22.5% 100.0% 99.3% 99.7% 92.6% 30.7% 99.8% 95.7% 60.6% 99.3% 3.9% 99.8% 7.8% 8.6% 99.9% 0.3% 49.6%

least with the reduction ratio of 0.3%. Additionally, most optimal values of Kset(NH4), Kd(NH4) and Rd(NH4) according to the new water quality series are very close to their initial optimal values. Furthermore, most of the optimal values of Rd(COD) and Rset(NH4) are close to the maximum values, most of the optimal values of Rset(COD) and Kd(COD) are close to the median values and those of Rset(orgN) are close to the minimum values. All of these optimal values are quite different from the initial optimal values due to the parameter uncertainties. Compared with the marginal final distributions due to uncertainty propagation of hydrological parameters, those of Kset(NH4), Kd(NH4), Rd(NH4) and Rset(orgN) are not changed obviously, which account for 50% of the total selected water quality parameters. However, the marginal final distributions of Rd(COD), Rset(COD) and Rset(NH4) are changed from the negative to positive skewness, while the change of Kd (COD) distribution is the opposite. 3.4.2. Uncertainty intervals of water quality simulation The water quality simulations due to the parameter uncertainties of water quality modules at the 95% confidence interval are shown in Fig. 8a. The uncertainty intervals are quite narrow with the ARIL of only 0.12. Compared with the observations, only some observations in the flood season (June–September) are bracketed in the uncertainty intervals, and account for 12.1% of total observations (i.e., P95CI = 12.1%) (Table 4). All the simulation groups at the 95% confidence interval are underestimated, particularly in the nonflood season (October–May) with the ranges of Nbias from 0.27 to 0.44. The temporal variations of simulations match well with the observations (i.e., the ranges of r are from 0.72 to 0.77), but the NS is slightly weakened with the ranges from 0.23 to 0.46. The uncertainties at the 95% confidence interval become larger when the parameter uncertainties of hydrological cycle modules are involved (Fig. 8b). The ARIL is 0.41, much greater than that only caused by the parameter uncertainties of water quality modules (i.e., ARIL = 0.12). The most obvious changes appear in the nonflood season, particular the pre-flood season from January to May. Furthermore, the P-95CI is also increased to 21.0%, particularly in January and February due to the wider uncertainty intervals. The model performance criteria

Fig. 4. Daily runoff observations, optimal simulations and uncertainty intervals at the 95% confidence level due to parameter uncertainty (a), parameter and module structure uncertainties (b).

negative skewness with the values from −0.175 to −1.278. Compared with the initial ranges of parameter values, all the ranges of final distributions at the 95% confidence interval are shortened remarkably. The changes of Kset(NH4), Kd(NH4) and Rd(NH4) are the greatest whose reduction ratios are close to 100%, while the change of Rset(NH4) is the 630

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Table 4 Assessments of uncertainty intervals of runoff and water quality simulations due to uncertainties of parameters, module structures and parameter uncertainty propagations. Variables

Runoff simulation

Water quality simulation due to uncertainty propagation of hydrological parameters Water quality simulation

Uncertainty sources

Parameters of hydrological cycle modules Module structure and parameters of hydrological cycle modules Parameters of hydrological cycle modules Parameters of hydrological cycle modules and module structure of water quality modules Parameters of water quality modules Parameters of hydrological cycle and water quality modules Parameters of hydrological cycle and water quality modules, module structure of water quality modules

Range of Nbias

Range of r

Range of NS

Other criteria

Min

Max

Min

Max

Min

Max

ARIL

P-95CI

−0.08 −0.10

0.01 0.04

0.89 0.84

0.90 0.89

0.79 0.70

0.81 0.78

0.37 3.66

18.7% 86.0%

0.17 −0.24

0.23 0.12

0.75 0.24

0.78 0.72

0.48 −0.43

0.55 0.46

0.05 5.94

1.7% 94.8%

0.27 0.27

0.44 0.68

0.72 0.59

0.77 0.78

0.23 −0.25

0.46 0.46

0.12 0.41

12.1% 21.0%

−0.09

0.46

0.07

0.64

−0.52

0.38

5.94

92.0%

Fig. 5. Marginal posterior densities of all the selected water quality parameters of HEQM due to the uncertainty propagation of hydrological parameters. 631

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

uncertainty analysis are able to be resampled efficiently without extensive programming by the Bootstrap package of the R platform and are further off-line coupled with HEQM. Furthermore, it does not need to determinate the subjective acceptance threshold of likelihood function compared with the GLUE, which is also a widely used approach (Blasone et al., 2008; Freni et al., 2008; Li et al., 2010a). Although this approach is conducted based on a great number of model iterations, the parallel optimization is advantageous to enhance the computation efficiency (Li et al., 2010b; Shao et al., 2012). The results present that 18.7% of the whole runoff simulation uncertainties are from the parameter uncertainties of hydrological cycle module of HEQM, and the figure increases to 86.0% if the uncertainties of module structure are considered further. Although a majority of uncertainties are probably from the model structures, the parameter estimation is also a critical uncertainty source, which causes the inaccurate simulation performance of runoff to some extent. Our finding is basically in agreement with many existing studies. For example, the uncertainties from parameter, parameter and model structure contribute 18.8% and 82.3% of the total runoff simulation uncertainties for DTVGM (Distributed Time Variant Gain Model), respectively in the Chaobai River basin (Wang et al., 2009; Li et al., 2010a), and 12.1–44.1% and 79.2% for WASMOD (Water And Snow balance MODelling system), respectively in the central Sweden (Engeland et al., 2005), Chao River basin (Li et al., 2010a) and Shuntian basin (Jin et al., 2010). Additionally, the parameter uncertainty contributes 10.0–13.0% of total runoff simulation uncertainties for SWAT in the Chao River basin (Yang et al., 2007) and Yingluoxia watershed (Li et al., 2010b). In consideration of the interactions between hydrological and water quality processes, the uncertainties of hydrological parameters are inevitable to cause the inaccurate simulations of subsequent water quality concentrations (Freni et al., 2008; Quinton et al., 2010). Our results further highlight this viewpoint by the uncertainty propagation analysis and show that the uncertainties of hydrological parameters introduce a certain simulation uncertainties of water quality by changing the final distributions of water quality parameters. The most affected parameters are the deoxygenation or settling rates of water quality variables (e.g., Rd(COD), Rset(orgN) and Kd(COD)), all of which are highly related with the runoff magnitude and velocity (Adrian and Sanders, 1998; Marttila and Klove, 2010; Freni et al., 2011; Jung et al., 2015). Although the module structure is the major uncertainty source, the uncertainty propagation of hydrological parameters still slightly weakens the simulation performance of water quality concentration. The uncertainties of water quality parameters and module structures themselves also contribute a proportion of simulation uncertainties of water quality. For SWAT, 71–85% of the total sediment simulation uncertainties at the 95% confidence interval (Abbaspour et al., 2007; Talebizadeh et al., 2010) and over 75% of total phosphorous simulation uncertainties at the 90% confidence interval (Gong et al., 2011) are caused by the uncertainties of parameters and model structures. Due to the ill-posed water quality modules of the stormwater model SIMPLE KAREN, only 22.8–28.8% of total simulation uncertainties of total suspended solids at the 90% confidence interval are from the uncertainties of parameters and model structures (Dotto et al., 2012). Our study investigates the uncertainties of parameters and module structures step-by-step and finds that 12.1% of total NH4-N simulation uncertainties are from the parameter uncertainties of water quality modules, and 92.0% are from the uncertainties of parameters and module structures which are very close to the optimal uncertainties (95%). Thus, the results calculated by our proposed framework are much more satisfying than those of above-mentioned studies. Additionally, the uncertainty propagation of module structure is not investigated in our study because besides the river flow routing, numerous of other hydrological processes and their corresponding parameters are able to affect the simulation performance of water quality, such as surface runoff generation, soil water and its flow, infiltration (Tesoriero et al., 2009; Fiener et al., 2011; Parajuli et al., 2013). In this

Fig. 6. Daily observations of water quality concentrations, optimal simulations and uncertainty intervals at the 95% confidence level due to parameter propagation uncertainty (a), parameter and module structure uncertainties (b).

range from 0.27 to 0.68 for Nbias, from 0.59 to 0.78 for r, and from −0.25 to 0.46 for NS. All the simulation groups at the 95% confidence interval are still underestimated and the temporal variations of simulations do not change obviously, but the minimum NS decreases a lot. Moreover, due to further consideration of structure uncertainties of water quality modules, the 95% confidence interval becomes much wider (Fig. 8c) and the ARIL further increases to 5.94. Moreover, most of total observations are bracketed, only except some high observations in April of the nonflood season. The P-95CI is 92.0%, which is also close to the optimum of P-95CI (95%). Thus, the uncertainty intervals of water quality simulations due to the uncertainties of parameters and module structures are very robust by our proposed integrated framework. Although both ARIL and P-95CI are much greater than those caused by parameter uncertainties of water quality modules (0.12 and 21.0%), and parameter uncertainties of both hydrological and water quality modules (0.41 and 92.0%), they are very close to or even the same with these values due to the uncertainties of parameters of hydrological cycle modules, and structure of water quality modules. Compared with the effect of uncertainties of module structures, the effect of parameter uncertainties is much smaller, particularly for the parameter uncertainty propagation of hydrological cycle module. However, the model performance becomes weaker with the ranges of Nbias from −0.09 to 0.46, the ranges of r from 0.07 to 0.64, and the ranges of NS from −0.52 to 0.38. 4. Discussion Our study proposes an assessment framework of uncertainty and its propagation for integrated water system models by coupling the frequentist approach with model auto-calibration techniques (i.e., the Bootstrap method and SCE-UA), both of which are robust and widely adopted in the uncertainty analysis and model calibration, respectively (Efron, 1979; Duan et al., 1992, 1994). It is much easier to implement than the Bayesian approach because the model residuals used for the 632

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Fig. 7. Marginal posterior densities of all the selected water quality parameters of HEQM due to the parameter uncertainty of water quality modules.

5. Conclusions

study, all the selected parameters of hydrological cycle module are about the soil water, runoff generation and baseflow in the sub-catchments, rather than the river routing parameters based on the parameter sensitivity analysis. Although the runoff simulation residuals are deduced, the residuals of the other hydrological variables (e.g., soil water, evapotranspiration) are not obtained due to the observation limitations. Thus, the uncertainty intervals of runoff simulation are only a part of the whole structure uncertainties of hydrological cycle modules, and the others are still not available for the assessment of uncertainty propagation effects on the simulation performance of water quality. In next work, the assessment of uncertainty contributions of specific hydrological or water quality processes should be conducted in details to support the investigations of the uncertainty propagation of module structure in the integrated models.

Multiple uncertainties of a complicated integrated water system are investigated by integrating the Bootstrap resampling method and autocalibration algorithm (SCE-UA), including uncertainties of parameter and structure of hydrological cycle modules, parameter uncertainty propagation to the subsequent water quality modules, as well as uncertainties of parameter and structure of water quality modules. Results show that: HEQM is able to simulate the daily runoff very well, and the NH4-N concentration satisfactorily by the auto-calibration. The uncertainty interval of runoff simulation due to parameter uncertainties is narrow and only 18.7% of all the observations fall into the uncertainty intervals at the 95% confidence level. The uncertainty intervals due to the uncertainties of parameter and module structure become much 633

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

mechanism investigation should be strengthened and the empirical or simplified equations should be replaced by the physically or mathematically based equations. For instance, more types of land use/land cover units should be designed in the hydrological cycle module of HEQM to relieve the impacts of spatial heterogeneities based on the high spatial resolutions of vegetation coverage. River sediment resuspension and deposition processes should be introduced in the water quality module to improve the simulation of water quality, particularly during the flood season in the regulated rivers. Acknowledgements This study was supported by Natural Science Foundation of China (No. 41671024), the Program for “Bingwei” Excellent Talents in Institute of Geographic Sciences and Natural Resources Research, CAS, China (No. 2015RC201), the International Fellowship Initiative, Institute of Geographic Sciences and Natural Resources Research, CAS, China (No. 2017VP04), and the China Youth Innovation Promotion Association CAS, China (No. 2014041). Thanks also to the editors, and three anonymous referees for their constructive comments. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.jhydrol.2018.08.070. References Abbaspour, K.C., Yang, J., Maximov, I., Siber, R., Bogner, K., Mieleitner, J., Zobrist, J., Srinivasan, R., 2007. Modelling hydrology and water quality in the pre-alpine/alpine Thur watershed using SWAT. J. Hydrol. 333 (2), 413–430. Adrian, D.D., Sanders, T.G., 1998. Oxygen sag equation for second-order bod decay. Water Res. 32 (3), 840–848. Arnold, J.G., Srinivasan, R., Muttiah, R.S., Williams, J.R., 1998. Large-area hydrologic modeling and assessment: Part I. Model development. J. Am. Water Resour. Assoc. 34, 73–89. Arsenault, R., Poissant, D., Brissette, F., 2015. Parameter dimensionality reduction of a conceptual model for streamflow prediction in Canadian, snowmelt dominated ungauged basins. Adv. Water Resour. 85, 27–44. Bates, B.C., Campbell, E.P., 2001. A Markov chain Monte Carlo scheme for parameter estimation and inference in conceptual rainfall–runoff modeling. Water Resour. Res. 37 (4), 937–947. Beck, M.B., 1987. Water quality modeling: a review of the analysis of uncertainty. Water Resour. Res. 23 (8), 1393–1442. Becker, R.A., Chambers, J.M., Wilks, A.R., 1988. The New S Language. Wadsworth & Brooks/Cole. Beven, K.J., Binley, A., 1992. Future of distributed models: model calibration and uncertainty prediction. Hydrol. Process. 6 (3), 279–298. Beven, K.J., Freer, J., 2001. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology. J. Hydrol. 249 (1–4), 11–29. Blasone, R.S., Madsen, H., Rosbjerg, D., 2008. Uncertainty assessment of integrated distributed hydrological models using GLUE with Markov chain Monte Carlo sampling. J. Hydrol. 353 (1–2), 18–32. China’s national standard Current land use condition classification (GB/T21010–2007), 2007. General administration of quality supervision, inspection and quarantine of China and Standardization administration of China, Beijing, China (In Chinese). Crosetto, M., Ruiz, J.A.M., Crippa, B., 2001. Uncertainty propagation in models driven by remotely sensed data. Remote Sens. Environ. 76 (3), 373–385. Dai, Y., Zeng, X., Dickinson, R.E., Baker, I., Bonan, G.B., Bosilovich, M.G., Denning, A.S., Dirmeyer, P.A., Houser, P.R., Niu, G., Oleson, K.W., Schlosser, C.A., Yang, Z.L., 2003. The common land model. Bull. Amer. Meteor. Soc. 84, 1013–1023. Dotto, C.B., Mannina, G., Kleidorfer, M., Vezzaro, L., Henrichs, M., Mccarthy, D.T., Freni, G., Rauch, W., Deletic, A., 2012. Comparison of different uncertainty techniques in urban stormwater quantity and quality modelling. Water Res. 46 (8), 2545–2558. Duan, Q., Sorooshian, S., Gupta, V.K., 1994. Optimal use of the SCEUA global optimization method for calibrating watershed models. J. Hydrol. 158, 265–284. Duan, Q.Y., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall–runoff models. Water Resour. Res. 28 (4), 1015–1031. Efron, B., 1979. Bootstrap methods: another look at the jackknife. Ann. Statist. 7, 1–26. Engeland, K., Gottschalk, L., 2002. Bayesian estimation of parameters in a regional hydrological model. Hydrol. Earth Syst. Sci 6 (5), 883–898. Engeland, K., Xu, C.Y., Gottschalk, L., 2005. Assessing uncertainties in a conceptual water balance model using Bayesian methodology. Hydrol. Sci. J. 50 (1), 45–63. Fiener, P., Auerswald, K., Oost, K.V., 2011. Spatio-temporal patterns in land use and management affecting surface runoff response of agricultural catchments—a review. Earth-Sci. Rev. 106 (1–2), 92–104.

Fig. 8. Daily observations of water quality concentrations, optimal simulations and uncertainty intervals at the 95% confidence level due to parameter uncertainties of water quality modules (a), parameter uncertainties of both hydrological cycle and water quality modules (b), uncertainties of parameters of both hydrological cycle and water quality modules, module structure of water quality modules (c).

wider and 86.0% of all the observations (mainly excluding some high flows) fall into the 95% confidence interval. Moreover, the uncertainty propagation of hydrological parameters does not obviously change the simulated variations of water quality concentrations throughout the period, and only 1.7% of observations fall into the uncertainty intervals. However, the simulation uncertainty intervals are further increased obviously due to the parameter uncertainties of water quality modules, and the proportions of observations falling into the confidence interval are increased from 12.1% to 21.0%. Along with the consideration of structure uncertainties of water quality modules, the intervals are much wider and contain 92.0% of total observations, only except the extreme high values. The module structure error is still a major uncertainty source for both runoff and water quality simulations, although the parameter uncertainty and its propagation cannot be ignored. The model structure improvement is still the priority for the integrated water system modelling and can greatly improve the simulation performance of individual modules and reduce their uncertainty propagation. More

634

Journal of Hydrology 565 (2018) 623–635

Y. Zhang, Q. Shao

Paola, C., Foufoula-Georgiou, E., Dietrich, W.E., Hondzo, M., Mohrig, D., Parker, G., Power, M.E., Rodriguez-Iturbe, I., Voller, V., Wilcock, P., 2006. Toward a unified science of the Earth’s surface: opportunities for synthesis among hydrology, geomorphology, geochemistry, and ecology. Water Resour. Res. 42, W03S10. https:// doi.org/10.1029/2005WR004336. Parajuli, P.B., Jayakody, P., Sassenrath, G.F., Ouyang, Y., Pote, J.W., 2013. Assessing the impacts of crop-rotation and tillage on crop yields and sediment yield using a modeling approach. Agric. Water Manage. 119 (119), 32–42. Quinton, J.N., Govers, G., Oost, K.V., Bardgett, R.D., 2010. The impact of agricultural soil erosion on biogeochemical cycling. Nat. Geosci. 3 (5), 311–314. R Development Core Team, 2010. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria http://www.Rproject.org. Refsgaard, J.C., Sluijs, J.P.V.D., Brown, J., Keur, P.V.D., 2006. A framework for dealing with uncertainty due to model structure error. Adv. Water Resour. 29 (11), 1586–1597. Shao, Q.X., Lerat, J., Brink, H., Tomkins, K., Yang, A., Peeters, L., Li, M., Zhang, L., Podger, G., Renzullo, L.J., 2012. Gauge based precipitation estimation and associated model and product uncertainties. J. Hydrol. 444–445, 100–112. Shao, Q.X., Lerat, J., Podger, G., Dutta, D., 2014. Uncertainty estimation with bias-correction for flow series based on rating curve. J. Hydrol. 510 (3), 137–152. Sharpley, A. N., Williams, J. R., 1990. EPIC-erosion/productivity impact calculator: 1. Model documentation. Technical Bulletin- United States Department of Agriculture, Agric. Res. Service, Washington D.C., USA. Talebizadeh, M., Morid, S., Ayyoubzadeh, S.A., Ghasemzadeh, M., 2010. Uncertainty analysis in sediment load modeling using ANN and SWAT model. Water Resour. Manage. 24 (9), 1747–1761. Tesoriero, A.J., Duff, J.H., Wolock, D.M., Spahr, N.E., Almendinger, J.E., 2009. Identifying pathways and processes affecting nitrate and orthophosphate inputs to streams in agricultural watersheds. J. Environ. Qual. 38, 1892–1900. Todini, E., 2007. Hydrological catchment modelling: past, present and future. Hydrol. Earth Syst. Sci. 11 (1), 468–482. Wang, G., Xia, J., Chen, J., 2009. Quantification of effects of climate variations and human activities on runoff by a monthly water balance model: a case study of the Chaobai River basin in northern China. Water Resour. Res. 45, W00A11. https://doi. org/10.1029/2007WR006768. Willems, P., 2008. Quantification and relative comparison of different types of uncertainties in sewer water quality modeling. Water Res. 42 (13), 3539–3551. Wucchini, W., Adamson, P.T., 1988. On the application of the Bootstrap to assess the risk of deficient annual inflows to a reservoir. Water Resour. Manag. 2 (4), 245–254. Xu, C.Y., Tunemar, L., Chen, Y.D., Singh, V.P., 2006. Evaluation of seasonal and spatial variations of conceptual hydrological model sensitivity to precipitation data errors. J. Hydrol. 324, 80–93. Yang, J., Reichert, P., Abbaspour, K.C., 2007. Bayesian uncertainty analysis in distributed hydrological modelling: a case study in the Thur River basin (Switzerland). Water Resour. Res. 43, W10401. https://doi.org/10.1029/2006WR005497. Yen, H., Wang, X., Fontane, D.G., Harmel, R.D., Arabi, M., 2014. A framework for propagation of uncertainty contributed by parameterization, input data, model structure, and calibration/validation data in watershed modeling. Environ. Modell. Softw. 54 (3), 211–221. Yue, S., Wang, C.Y., 2002. Assessment of the significance of sample serial correlation by the bootstrap test. Water Resour. Manag. 16 (1), 23–35. Zhang, Y.Y., Shao, Q.X., Taylor, J.A., 2016a. A balanced calibration of water quantity and quality by multi-objective optimization for integrated water system model. J. Hydrol. 538 (538), 802–816. Zhang, Y.Y., Shao, Q.X., Ye, A.Z., Xing, H.T., Xia, J., 2016b. Integrated water system simulation by considering hydrological and biogeochemical processes: model development, with parameter sensitivity and autocalibration. Hydrol. Earth Syst. Sci. 20, 529–553. Zhang, Y.Y., She, D.X., Xia, J., Zhai, X.Y., Zuo, Q.T., 2015. Causal analysis on the specified paroxysmal water pollution incidents in Huai River Basin. Environ. Eng. Manag. J. 14 (1), 139–151. Zucchini, W., Adamson, P.T., 1989. Bootstrap confidence intervals for design storms from exceedance series. Hydrol. Sci. J. 34 (1), 41–48.

Fowler, K., Peel, M., Western, A., Zhang, L., Peterson, T., 2016. Simulating runoff under changing climatic conditions: revisiting an apparent deficiency of conceptual rainfallrunoff models. Water Resour. Res. 52 (3), 1820–1846. Freni, G., Mannina, G., Viviani, G., 2011. Assessment of the integrated urban water quality model complexity through identifiability analysis. Water Res. 45 (1), 37–50. Freni, G., Mannina, G., Viviani, G., 2008. Uncertainty in urban stormwater quality modelling: the effect of acceptability threshold in the GLUE methodology. Water Res. 42 (8–9), 2061–2072. Freni, G., Mannina, G., Viviani, G., 2009. Assessment of data availability influence on integrated urban drainage modelling uncertainty. Environ. Modell. Softw. 24 (10), 1171–1181. Gabellani, S., Boni, G., Ferraris, L., Hardenberg, J.V., Provenzale, A., 2007. Propagation of uncertainty from rainfall to runoff: a case study with a stochastic rainfall generator. Adv. Water Resour. 30 (10), 2061–2071. Gallagher, M., Doherty, J., 2007. Parameter estimation and uncertainty analysis for a watershed model. Environ. Modell. Softw. 22 (7), 1000–1020. Gong, Y., Shen, Z., Hong, Q., Liu, R., Liao, Q., 2011. Parameter uncertainty analysis in watershed total phosphorus modeling using the GLUE methodology. Agric. Ecosyst. Environ. 142 (3), 246–255. GWSP, 2005. The Global Water System Project: Science Framework and Implementation Activities. Earth System Science Partnership. Hamrick, J.M., 1992. A Three-Dimensional Environmental Fluid Dynamics Computer Code: Theoretical and Computational Aspects, Special Report, The College of William and Mary. Virginia Institute of Marine Science, Virginia, USA, pp. 317. Harmel, R.D., King, K.W., Haggard, B.E., Wren, D.G., Sheridan, J.M., 2006. Practical guidance for discharge and water quality data collection on small watersheds. Trans. ASABE 49 (4), 937–948. Hirsch, R.M., Archfield, S.A., Cicco, L.A.D., 2015. A bootstrap method for estimating uncertainty of water quality trends. Environ. Modell. Softw. 73 (C), 148–166. Hu, Y.M., Liang, Z.M., Liu, Y.W., Zeng, X.F., Wang, D., 2015. Uncertainty assessment of estimation of hydrological design values. Stoch. Env. Res. Risk A 29 (2), 501–511. Jia, H., Xu, T., Liang, S., Zhao, P., Xu, C., 2018. Bayesian framework of parameter sensitivity, uncertainty, and identifiability analysis in complex water quality models. Environ. Modell. Softw. 104, 13–26. Jin, X., Xu, C.Y., Zhang, Q., Singh, V.P., 2010. Parameter and modeling uncertainty simulated by GLUE and a formal Bayesian method for a conceptual hydrological model. J. Hydrol. 383 (3), 147–155. Jung, B.J., Jeanneau, L., Alewell, C., Kim, B., Park, J.H., 2015. Downstream alteration of the composition and biodegradability of particulate organic carbon in a mountainous, mixed land-use watershed. Biogeochemistry 122 (1), 79–99. Leta, O.T., Nossent, J., Velez, C., Shrestha, N.K., Griensven, A.V., Bauwens, W., 2015. Assessment of the different sources of uncertainty in a swat model of the river senne (Belgium). Environ. Modell. Softw. 68 (C), 129–146. Li, C., Frolking, S., Frolking, T.A., 1992. A model of nitrous oxide evolution from soil driven by rainfall events: 1. Model structure and sensitivity. J. Geophys. Res. 97 (D9), 9759–9776 (1984–2012). Li, L., Xia, J., Xu, C.Y., Singh, V.P., 2010a. Evaluation of the subjective factors of the GLUE method and comparison with the formal Bayesian method in uncertainty assessment of hydrological models. J. Hydrol. 390 (3–4), 210–221. Li, Z., Shao, Q., Xu, Z., Cai, X., 2010b. Analysis of parameter uncertainty in semidistributed hydrological models using bootstrap method: a case study of SWAT model applied to Yingluoxia watershed in northwest China. J. Hydrol. 385 (1–4), 76–83. Liang, X., Lettenmaier, D.P., Wood, E.F., Burges, S.J., 1994. A simple hydrologically based model of land surface water and energy fluxes for general circulation models. J. Geophys. Res. 99 (D7), 14415–14428. Marttila, H., Klove, B., 2010. Managing runoff, water quality and erosion in peatland forestry by peak runoff control. Ecol. Eng. 36 (7), 900–911. Montanari, A., Brath, A., 2004. A stochastic approach for assessing the uncertainty of rainfall–runoff simulations. Water Resour. Res. 40, W01106. Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Binger, R.L., Harmel, R.D., Veith, T., 2007. Model evaluation guidelines for systematic quantification of accuracy in watershed simulations. Trans. ASABE 50 (3), 885–900. Novic, A.J., Ort, C., O'Brien, D.S., Lewis, S.E., Davis, A.M., Mueller, J.F., 2018. Understanding the uncertainty of estimating herbicide and nutrient mass loads in a flood event with guidance on estimator selection. Water Res. 132, 99–110.

635