Addressing ten questions about conceptual rainfall–runoff models with global sensitivity analyses in R

Journal of Hydrology 503 (2013) 135–152

Contents lists available at ScienceDirect

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

Addressing ten questions about conceptual rainfall–runoff models with global sensitivity analyses in R Mun-Ju Shin a,c,1, Joseph H.A. Guillaume b,c,1, Barry F.W. Croke a,b,c,2, Anthony J. Jakeman b,c,⇑ a

Department of Mathematics, Building 48A, Linnaeus Way, The Australian National University, Canberra, ACT 0200, Australia Integrated Catchment Assessment and Management Centre, Fenner School of Environment and Society, Building 48A, Linnaeus Way, The Australian National University, Canberra, ACT 0200, Australia c National Centre for Groundwater Research and Training, Fenner School of Environment and Society, Building 48A, Linnaeus Way, The Australian National University, Canberra, ACT 0200, Australia b

a r t i c l e

i n f o

Article history: Received 22 February 2013 Received in revised form 2 August 2013 Accepted 23 August 2013 Available online 7 September 2013 This manuscript was handled by Geoff Syme, Editor-in-Chief, with the assistance of Bellie Sivakumar, Associate Editor Keywords: Sensitivity analysis Rainfall–runoff model Identifiability

s u m m a r y Sensitivity analysis (SA) is generally recognized as a worthwhile step to diagnose and remedy difficulties in identifying model parameters, and indeed in discriminating between model structures. An analysis of papers in three journals indicates that SA is a standard omission in hydrological modeling exercises. We provide some answers to ten reasonably generic questions using the Morris and Sobol SA methods, including to what extent sensitivities are dependent on parameter ranges selected, length of data period, catchment response type, model structures assumed and climatic forcing. Results presented demonstrate the sensitivity of four target functions to parameter variations of four rainfall–runoff models of varying complexity (4–13 parameters). Daily rainfall, streamflow and pan evaporation data are used from four 10-year data sets and from five catchments in the Australian Capital Territory (ACT) region. Similar results are obtained using the Morris and Sobol methods. It is shown how modelers can easily identify parameters that are insensitive, and how they might improve identifiability. Using a more complex objective function, however, may not result in all parameters becoming sensitive. Crucially, the results of the SA can be influenced by the parameter ranges selected. The length of data period required to characterize the sensitivities assuredly is a minimum of five years. The results confirm that only the simpler models have well-identified parameters, but parameter sensitivities vary between catchments. Answering these ten questions in other case studies is relatively easy using freely available software with the Hydromad and Sensitivity packages in R. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction All too often the reporting of model evaluation exercises omits a key step, sensitivity analysis (SA), in helping to understand the performance of that model. Sensitivity analysis is the study of how the variation in the output of a model can be qualitatively or quantitatively apportioned to different sources of variation in the model input data, including model parameter values (Saltelli et al., 2000). It can be one of the simplest aids in diagnosing and remedying poor identifiability of models, to allow parameters to be more reliably estimated. However, as we illustrate below, when such analysis ⇑ Corresponding author at: Integrated Catchment Assessment and Management Centre, Fenner School of Environment and Society, Building 48A, Linnaeus Way, The Australian National University, Canberra, ACT 0200, Australia. Tel.: +61 (0)2 6125 4742; fax: +61 (0)2 6125 8395. E-mail addresses: [email protected] (M.-J. Shin), joseph.guillaume@ anu.edu.au (J.H.A. Guillaume), [email protected] (B.F.W. Croke), tony. [email protected] (A.J. Jakema. 1 Tel.: +61 (0)2 6125 4670. 2 Tel.: +61 (0)2 6125 0666. 0022-1694/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2013.08.047

has been undertaken it is likely to have been a simple One-At-aTime (OAT) sensitivity to parameters over some local range. Saltelli and Annoni (2010) have previously shown local/single point OAT analysis to be a perfunctory and inadequate method of sensitivity analysis unless the model is linear. If the model parameters are nonlinear, local or single point OAT cannot sample the space of parameters adequately. Hence they suggest using global SA methods, such as the Morris (1991) or the Sobol (1990) variance-based method. Because the Sobol method has a larger computational cost, they recommend the Morris method when model runs are constrained by computational time. It is our impression that there has been a pervasive omission of SA in environmental model evaluation exercises, and that when SA has been invoked there has been an over-reliance on a local OAT approach. This is borne out by an illustrative analysis of rainfall– runoff modeling papers published in three journals in 2011: Journal of Hydrology, Ecological Modelling, and Environmental Modelling and Software. In 2011 these journals published 164 papers related to rainfall–runoff modeling (found when searching with the keywords: ‘‘rainfall–runoff’’ and ‘‘hydrological model’’ using

136

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

Science Direct, http://www.sciencedirect.com). Of the 164 papers in that year, there are only 59 that include the words ‘‘sensitivity analysis’’ and only 11 of those papers (7% of 164 papers) that use SA methods for their rainfall–runoff modeling. Of the 11 papers, two papers use a variance-based method, one paper uses the Morris screening method, five papers use an OAT method and the other three papers use a heuristic, Bayesian or multi-parameter nonlinear SA method that employs neural networks as a pseudo simulator to reduce the computational burden of the analysis. The remaining 48 papers report either sensitivity analysis using a percentage change of input data to quantify the impact of climate change or just mention that they either did a sensitivity analysis or did not; that is they either considered scenario-type changes in inputs (impact analysis rather than model evaluation), or did not present any results on application of SA at all. If OAT methods are removed from our count of papers, there are only 6 papers (4%) out of 164 papers that seem to use an adequate SA method to evaluate their rainfall–runoff models. It is clear therefore that model evaluation practices, as recently as 2011 for rainfall–runoff modeling, are demonstrably poor because even simple approaches of sensitivity analysis are not applied as a standard to investigate identifability. It is quite likely that such generally inadequate practices extend beyond rainfall–runoff modeling, not just in hydrology but in the environmental sciences more widely. Of course there are exceptions to our general contention that sensitivity analyses of models reported in the literature tend to be perfunctory. For example, van Griensven et al. (2006) pointed out that over-parameterisation is a well-known problem in rainfall–runoff models with large numbers of parameters. They showed that adopting the Morris method with Latin Hypercube sampling can reduce the number of parameters in the Soil and Water Assessment Tool (SWAT) model for more efficient use. They also concluded that there are clearly different SA results between the catchments, and therefore each catchment model requires its own sensitivity analysis to select a parameter set to be used for model calibration or uncertainty analysis. To indicate another exception, van Werkhoven et al. (2009) undertook a Sobol-based sensitivity analysis of the Sacramento model in order to identify and fix insensitive parameters that have a total sensitivity index less than some threshold. While reducing the complexity of the calibration process via SA, the resultant model had essentially the same predictive performance to the non-fixed parameter case. In this paper we show that an informative sensitivity analysis of a model need not require more than a few lines of code. Two SA methods, Morris (1991) and Sobol (1990), are used on four wellknown rainfall–runoff models with four target functions, and four decadal daily datasets from five study catchments for illustrative purposes. For the analysis the SA package ‘‘Sensitivity’’ in R is used. It includes some global SA methods including the Morris and Sobol algorithms applied in this study. This package is linked to the ‘‘Hydrological Model Assessment and Development’’ (Hydromad) package (Andrews et al., 2011), an open source software package which is also built on the R platform. Hydromad and its user guide are available from http://hydromad.catchment.org. It includes the rainfall–runoff models investigated here except for the SIMHYD model which is called from Hydromad and also allows one to insert code for any conceptual hydrological model so that it can be evaluated. The paper aims to show how sensitivity analysis can be easily used to help answer 10 important questions. All these questions are to some extent case-specific. While the method used is transferrable across case-studies, caution is usually required in transferring the results. In Section 2, we describe the two SA methods examined in this paper. Section 3 summarizes the catchments, input data, target functions and methods used for our analysis. In

Section 4, we provide some answers to the 10 questions using the two SA methods. Section 5 presents a discussion of the study conclusions and our suggestions for future work. In Appendix A, we present a description of each of the four rainfall–runoff models used. Appendix B summarizes the hydro-climatic characteristics of the five catchments used in this study. In Appendix C, we present the simple code for implementing the two SA methods in R. 2. Algorithms for sensitivity analysis Brief descriptions of the two SA algorithms are given below. More detailed information can be found in the references given. 2.1. The Morris screening method Morris (1991) developed a global screening OAT design to determine which parameters are (a) negligible, (b) linear and additive, (c) nonlinear or involved in interactions with other parameters. This is a global sensitivity method because it searches the entire range of the parameters (Cariboni et al., 2007; Saltelli et al., 2000). It produces the elementary effect or difference (di) as a local sensitivity measure. But because it averages the nominated number of local/partial effects computed at different points using evenly spaced points of a parameter in the entire parameter space, the final measure is a global measure (Norton, 2009). In changing the values of parameter values (x) OAT, the algorithm calculates target functions for each parameter set. The difference of target function values is divided by the difference of parameter values to obtain the di function of parameters (x) as follows:

di ðxÞ ¼

½yðx1 ; x2 ; . . . ; xi1 ; xi þ D; xiþ1 ; . . . ; xk Þ  yðxÞ D

ð1Þ

where D is a multiple of 1/(p  1) selected to sample p equally spaced parameter values, spanning the range of the parameter, and y(x) is target function value for the parameter values x. The di calculation process is repeated a number of times (r), and the mean (l) and standard deviation (r) values of the r samples of di are used as sensitivity indices. The l and r indicate the influence of each parameter on the target function. A high l value implies that a parameter has an important overall influence on the target function, and a high value of r implies that a parameter has strong interactions with other parameters (i.e. one parameter interacts with other parameters) or the effect of the parameter is nonlinear (Morris, 1991; van Griensven et al., 2006). Instead of using l, this study uses the mean of the absolute values of the r samples of di denoted as l as suggested by Campolongo et al. (2007). The l solves the problem of the effects of opposite signs due to a model’s non-monotonic characteristics. The advantages of this method are that it has relatively low computational cost, is simple to understand and implement, and easy to interpret its results (Campolongo and Saltelli, 1997; Saltelli et al., 2000). It is useful for sensitivity analyses that involve models with a large number of parameters (e.g. Cariboni et al., 2007). The disadvantage of the Morris method is that it cannot estimate individual interactions between parameters, thereby giving only the overall interaction of a parameter with the rest of the model (Saltelli et al., 2000). It produces relative sensitivity measures by ranking the input parameters in order of significance but cannot quantify in absolute terms how much one parameter is more important than another (Campolongo and Saltelli, 1997; Alam et al., 2004). Therefore this method is valuable for preliminary parameter sensitivity analyses to get an overview of parameter sensitivity (Campolongo and Saltelli, 1997). For more descriptions of the Morris algorithm, see Morris (1991) and Saltelli et al. (2000). Various studies using the Morris method include Campolongo and Saltelli (1997), van

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

Griensven et al. (2006), Campolongo et al. (2007), Cariboni et al. (2007), Yang (2011), and Sun et al. (2012). 2.2. The variance-based method of Sobol The Sobol method is a variance-based, quantitative sensitivity analysis algorithm published first in Russian (Sobol, 1990) and later translated into English (Sobol, 1993). This method decomposes the model output variance into relative contributions from individual parameters and parameter interactions, and as a result the sensitivity of a given parameter is quantified by the ratio (ranging from 0 to 1) of its contribution to the output variance (e.g. van Werkhoven et al., 2009). It generates a First-Order Sensitivity Index (FSI) and a Total Sensitivity Index (TSI). The FSI is defined as the partial contribution of each parameter to the output variance divided by the total output variance. It measures the main effect of each parameter on the output (Saltelli et al., 2000). It can be defined as (Saltelli and Annoni, 2010):

V X ðEX ðYjX i ÞÞ Si ¼ i i VðYÞ

ð2Þ

where Xi is the ith parameter and Xi means the vector of all parameters not including Xi The inner expectation operator denotes that the mean of Y, the scalar objective function values, is considered over all possible values of Xi while keeping Xi fixed. The outer variance in the numerator is considered over all possible values of Xi, and the variance V(Y) in the denominator is the total (unconditioned) variance. The numerator of Si can be interpreted as the expected reduction in variance that would be obtained with fixed Xi (Saltelli and Annoni, 2010). The TSI can be defined as:

STi ¼

EXi ðV Xi ðYjXi ÞÞ VðYÞ

ð3Þ

and the interpretation of the numerator of STi is that it is the expected variance that would be left if all parameters are fixed except for Xi (Saltelli and Annoni, 2010). As proven by Sobol (2001) and Alısß and Rabitz (2001), the variance decomposition of the Sobol method provides information on both the parameter contribution to the output variance and model structure (Borgonovo, 2006) but this method is unique only if the input variables are independent (Borgonovo, 2006; Kucherenko et al., 2012). Even though the variance decomposition of the Sobol method is not unique in the case of a nonlinear model with non-independent parameters, still the interactions of a parameter with others are taken into account and can be estimated relative to the total sensitivity of that parameter. This paper uses Saltelli’s scheme (Saltelli, 2002) which computes FSI and TSI at a reduced cost of n(k + 2) from n(2k + 2) with a Monte Carlo random sampling method for the input variables. Here, n represents the initial sample size used (10,000 in this study) for generating the two sets of Monte Carlo estimates, and k is the number of parameters of the respective hydrological model (see Section 3). The Sobol method generates new parameter sets by replacing each parameter in one sample set by the same parameter in the other sample set, one by one. Hence the total number of samples generated is 60,000 for four parameter models, and 110,000 and 150,000 for nine and 13 parameter models, respectively. The TSI yields the total effect of a parameter including all its interactions with other parameters, while the FSI yields the effect of parameter Xi by itself (Saltelli and Annoni, 2010). The TSI does not provide a perfect characterization of the system but gives more reliable results than the FSI in terms of investigation of the overall effect of each parameter on the output (Saltelli et al., 2000). The Sobol algorithm has high computational cost. It is therefore useful for

137

detailed sensitivity analysis of the remaining parameters after performing a preliminary sensitivity analysis using a screening method such as that of Morris (e.g. Cariboni et al., 2007). For more description of the Sobol algorithm, see Saltelli et al. (2000). There are various reported studies using the Sobol method (Homma and Saltelli, 1996; Yatheendradas et al., 2008; van Werkhoven et al., 2009; Nossent et al., 2011; Yang, 2011). The Morris l has strong correlation with the Sobol TSI, and Morris r with the difference between TSI and FSI (Campolongo and Saltelli, 1997). There could be a discrepancy between identification of parameters with both algorithms; for example several parameters that are well identified (sensitive) by Morris may not be well identified by Sobol and vice versa (Campolongo and Saltelli, 1997). Both methods may fail to identify a significant parameter; for example a significant parameter is concluded to be an insignificant parameter, but according to Campolongo and Saltelli (1997) there are no cases reported for the reverse problem. 3. Models, catchments, input data, target functions and methods used We selected four routinely-used lumped conceptual hydrological models for sensitivity analysis: IHACRES, GR4J, Sacramento and SIMHYD models. These are of varying complexity, from four to 13 parameters. A description of each of the rainfall–runoff models and their parameters are given in Appendix A. Five catchments in the Australian Capital Territory (ACT) region in Australia were selected for the SA. They represent different hydro-climatologic conditions, varying substantially from east to west according to rainfall amounts, rainfall seasonality and hydrological response. A summary of the hydro-climatologic conditions of the five catchments is given in Appendix B. Daily data from four decadal periods (the 1970s, 1980s, 1990s, and 2000s) were used for the sensitivity analysis. The choice of a 10-year period of data is motivated by the suggestions of Yapo et al. (1996), Anctil et al. (2004) and Kim et al. (2011). Yapo et al. (1996) remarked that approximately eight years of data are required for the Sacramento model, the most highly parameterized of our four models used in this paper, to obtain calibrations that are relatively insensitive to the period selected. Kim et al. (2011) found that parameters for the simpler IHACRES model reached consistently estimated values at around eight years. Anctil et al. (2004) showed that five years are generally sufficient with the equally simple GR4J model structure. It should be noted that the length of the calibration period will depend on the amount of information contained in the associated time series data (driven largely by the number, range, timing and frequency of events), and the complexity of the catchment response. We therefore also assess how sensitivities vary across four decades with inherently different climate forcing. This study calculates the sensitivity of four target functions to changes in parameter values. Target functions can either be objective functions, which are measures of performance calculated by comparison of modeled and observed streamflow, or prediction functions, which are scalar values calculated from the modeled streamflow time series that are provided to the end-user for their practical use. For each decade the three objective functions invoked are the modified Nash Sutcliffe Efficiency (NSE) proposed by Mathevet et al. (2006), NSElog, and 0.5NSE + 0.5NSElog. The prediction function used is the frequency of low flows. The modified Nash–Sutcliffe Efficiency (NSE) is defined as follows:

( NSE ¼

),( ) Pn Pn 2 ðQ obs;i  Q sim;i Þ2 i¼1 ðQ obs;i  Q sim;i Þ 1 þ 1  Pi¼1 Pn 2 2 n i¼1 ðQ obs;i  Q obs Þ i¼1 ðQ obs;i  Q obs Þ ð4Þ

138

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

where n is the number of time steps, Qobs,i is observed flow at time step i (daily), Q obs is the mean of observed flow, Qsim,i is simulated flow, and the sum is over all time steps, up to a 10-year period. NSE has the advantage of reducing the range of values to the set [1,1], rather [1,1] as produced by the Nash–Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970). This reduces the influence of large negative values without otherwise changing the interpretation of the objective function. The log transformed NSE, denoted as NSElog indicates that all the variable terms in Eq. (4) are log transformed, while 0.5NSE + 0.5NSElog is a weighted average of the first two objective functions. These three are used as illustrative objective functions for our SA. The prediction function is intended to be meaningful to the end-user of model predictions. In this case, the prediction function selected is intended to measure the frequency of low flows, during which, for instance, pumping should not occur as it would deprive the river ecosystem of necessary water. This function is calculated as the proportion of days when modeled flow is below the historically observed 20th percentile flow level, denoted as F20. The type of sensitivity analysis applied in the paper can be undertaken very easily using the Sensitivity and Hydromad packages in R. After installing the packages, the requirements are: defining the target function(s), selecting the hydrological model, linking to the data sets, and calling one or more of the SA methods. The code for implementing the two SA methods in R is provided in Appendix C. Beyond the raw output of the SA, we present results obtained by summarizing the TSI values over multiple periods using Box plots, and combining multiple models and catchments in a single plot by representing the minimum and maximum TSI values. We also present the elementary effect of the Morris method for comparison with Sobol TSI values. 4. Results We first provide some introductory SA results in terms of TSI values from applying the Sobol method. These will serve as background for answering our 10 questions in the following subsections. Figs. 1–4 aim to summarize a large number of results. The four figures correspond to the four different rainfall–runoff models, respectively IHACRES, GR4J, Sacramento and SIMHYD. Each figure in turn contains a number of plots for different target functions. Each of these plots shows points for each parameter of the model, for five catchments. The coordinates of the points for each catchment are determined by the maximum (horizontal axis) and minimum (vertical axis) of the TSI values calculated by running the sensitivity analysis separately on each of four decades of data. Each point is therefore a summary of the TSI for one parameter of one model structure across four decades in one catchment with one target function. Parameter numbers are used for the different models. To identify the names of parameters in the plot, refer to the corresponding parameter numbers in Table A1 in Appendix A. If the aim of the sensitivity analysis is to select non-influential parameters with respect to the target function(s) and perhaps to fix their values, then the Sobol TSI has been suggested as a reasonable measure to use. The task then is to select a threshold below which parameters could be regarded as insensitive. A threshold of 0.2 is adopted here as illustration, meaning that the parameter produces at least 20% of the variance of the target function value. Although the study of van Werkhoven et al. (2009) adopts this and higher values, this threshold is still somewhat arbitrary, so care needs to be taken especially when the TSI values of parameters are near the threshold. 4.1. Question 1: Is the result of a sensitivity analysis affected by the parameter ranges selected? Global sensitivity analysis techniques like those of Morris and Sobol require the selection of a range of parameter values to be

explored. Reducing or expanding the parameter range could affect the sensitivity index, and could lead to insensitive parameters becoming sensitive or vice versa. Wang et al. (2013) for example showed that different parameter ranges for the WOFOST crop growth model yields differences in sensitive parameters when using the EFAST sensitivity analysis method. We investigate this further by repeating some sensitivity analyses with multiple parameter ranges. One simpler model (IHACRES here) and one more complex model (SIMHYD) are used for Sobol sensitivity analysis using one of the wetter/western (Gingera) and one of the drier/eastern (Tinderry) catchments with the NSE objective function. The range of the f and ss parameters of the IHACRES model, and the range of the InfiShp and InflCoef parameters of the SIMHYD model are changed by ±50%. The results are shown in Table 1 in terms of TSI values. Note that the original range of the parameter f (0.5–1.3) is selected based on the study of Croke and Jakeman (2004), and the range of the parameter ss (10–500) is set after considering the size of the catchments used in this study hence their likely response time. The original ranges of the parameter InfiShp (0–10) and InflCoef (0–400) are taken from Podger (2004) and more details are described in Appendix A. The reduction or expansion of the upper-bound of parameter range by ±50% in this study is arbitrary for illustrative purpose. For the IHACRES model, if the f parameter range is decreased by 50% (comparing IHACRES Original and f: 0.5–0.9 in Table 1) then the TSI value of the f parameter decreases by up to 82% in the Tinderry and Gingera catchments over the four decadal periods between 1970 and 2010. Indeed the f parameter has a lower TSI value in all cases and this is less than the 0.2 threshold except for the 2000s period. On the other hand sq and vs become more sensitive, having their TSI values increased up to 106% and 89% respectively in the Tinderry and Gingera catchments over the four decadal periods. Conversely, if the f parameter range is increased by 50% (comparing IHACRES Original and f: 0.5–1.7 in Table 1) then the TSI values of the f parameter increase by up to 88%, while the TSI values of the sq and vs parameters decrease by up to 40% and 39% respectively. But the sq and vs parameters are still sensitive parameters having TSI values over the threshold value 0.2 except for the 2000s. Therefore, the analysis confirms that in the case of the IHACRES model TSI values of parameters can change substantially if the range of a parameter is changed. Consider now the result of changing the range of the insensitive ss parameter of the IHACRES model and its effect on the other parameters. The range of the ss parameter is decreased (comparing IHACRES Original and ss: 10–250 in Table 1) and increased (comparing IHACRES Original and ss: 10–1000 in Table 1) by 50% but the ss parameter remains insensitive. Moreover, the effect on the TSI values of the other parameters is minor except for the vs parameter in the case of Tinderry and the 2000s period of analysis (0.21, +25% without rounding) when increasing the range of the ss parameter. Note that the 2000s constitute a very dry decade as explained in Appendix B, hence the sensitivity performance using data from the 2000s is different than in the other periods; when the range of the insensitive parameter is changed then it has only a minor effect on the TSI values of the parameter itself and on the other parameters. Similar results can be found with the more complex SIMHYD model. When the sensitive InfilShp parameter range is changed by ±50% (comparing SIMHYD Original, InfilShp: 0–5 and InfilShp: 0–20 in Table 1) then there are some changes (up to 32%) in the TSI value of the sensitive PercFrac parameter and significant changes (up to 521%) in the InfilShp parameter over the Gingera and Tinderry catchments. However this change of parameter range imparts only minor effect on the other insensitive six parameters and one somewhat sensitive InflCoef parameter, except for the case

139

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

0.5

(a)

1.0

4

3

(b) 1

Catchment Brindabella Burbong

Min. TSI

0.3

3 4

1 1

1

4 0.2

1

0.8

1

Gingera Orroral Crossing Tinderry

Min. TSI

0.4

1

0.6

1

0.4

4 4 3 0.1

0.0

3

2 2 22 2 0.0

1 3

0.0 0.2

3 4

0.2

0.4

0.6

0.8

1.0

3

3 4 3 22 3 22 2 0.0

4 4

0.2

0.4

Max. TSI 1.0

0.10

1

1

0.6

0.0

(d) 4

0.06

2

0.2

2 2

4

3 3

0.00

0.4

2

3

0.04

0.02

0.0

1.0

2

1

0.4

4 2 4 3 4 2 3 2 23 2 33

0.8

0.08

1 1 Min. TSI

Min. TSI

0.8

4

0.6

Max. TSI

(c)

0.2

1

4

0.6

0.8

3

0.00

1.0

0.05

0.10

0.15

0.20

Max. TSI

Max. TSI 0.8

1

(e)

0.10

(f)

1 1

0.08

2

1 Min. TSI

Min. TSI

0.6

0.4

4

1

0.06

2

0.04

2

0.2

0.0

2

4 4

2 2 2 2 4 4 33 2 3 33 0.0

0.2

0.02

0.00

0.4

0.6

Max. TSI

0.8

1.0

3 3 3 3 0.00

3 0.05

0.10

0.15

0.20

Max. TSI

Fig. 1. Plot of Sobol Total Sensitivity Indices (TSI) for the different target functions and IHACRES model parameters for the five catchments. Minimum and maximum TSI are calculated from the four decadal values. Note (a) is for NSE, (b) 0.5NSE + 0.5NSElog, (c) NSElog, (d) zoom of (c), (e) F20, and (f) zoom of (e). One-to-one line is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

140

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152 0.5

(a)

3

Catchment

2

Brindabella Burbong

Min. TSI

0.3

11

2

0.5

2

1

2

Gingera Orroral Crossing Tinderry

(b)

1

0.4

3

3 2

3 0.2

3 2

2 2 3

2 3

Min. TSI

0.4

1

3

2 1

3

0.3

3 0.2

1

1

1

0.1

4 4 0.0

4 4 0.0

0.0

0.2

1

0.1

4

0.4

0.6

0.8

Max. TSI

44 4 4 0.0

4 0.2

0.4

0.6

0.8

Max. TSI

0.6

(c) 0.5

3

2 2

2

3

Min. TSI

0.4

3 2

0.3

0.2

1 1

0.1

1

3 1

3 1 2

0.0

44 4 4 4 0.0

0.2

0.4

0.6

0.8

1.0

Max. TSI Fig. 2. Plot of Sobol Total Sensitivity Indices (TSI) for the different target functions and GR4J parameters for the five catchments. Minimum and maximum TSI are calculated from the four decadal values. Note (a) is for NSE, (b) NSElog, and (c) F20. One-to-one line is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of InflCoef parameter for Gingera with 50% increment of InfilShp parameter range. However in contrast to IHACRES, changing the range of sensitive parameters in the SIMHYD model has different effects on some parameter sensitivities for the wet versus the dry catchment datasets. The increased range of the InfilShp parameter (comparing SIMHYD Original and InfilShp: 0–20) increases the TSI value of the PervFrac parameter for Gingera (wet catchment) but decreases the TSI value of the PervFrac parameter for Tinderry (dry catchment), except for the 2000s period. The decreased range of the Inflcoef parameter (comparing SIMHYD Original and Inflcoef: 0–200) increases the TSI value of the PervFrac parameter for Gingera but decreases the TSI value of the PervFrac parameter for Tinderry (except for 2000s). Therefore different catchments and their response characteristics can have significant impact on sensitivity analysis of the SIMHYD model. This could well be caused by an interdependence of parameters. If this interdependence between parameters is ignored, then the models identified may perform poorly in different conditions. This supports the need to use SA methods that consider the interactions between parameters.

As the sensitivity of parameters can be strongly influenced by the range of parameter values used, it is important that the ranges used yield parameter sets that are considered plausible, such that results are not biased by implausible model realisations. Knowledge of the uncertainty in the parameter values is therefore needed to determine which model realisations to include. As this is not currently common practice, more work needs to be undertaken to assist in determining parameter ranges. The investigation of the response surface shapes of parameters due to the non-stationarity of parameter sensitivity could be a complementary study. As an interim measure, comparing sensitivities from multiple parameter ranges, as we have done here, may help determine whether results are robust to parameter range uncertainty by beginning to characterize the uncertainty in the sensitivity. 4.2. Question 2: Is the selected target function(s) insensitive to changes in any parameters? Table 2 shows the insensitive parameters in each model for three different target functions, NSE, NSElog and F20, for simple

141

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152 0.4

(a)

0.10

(b)

Catchment Brindabella Burbong

8 8

5

0.2

2

2

0.0

4

1 4 12 2 1 1 2 12111213 9912 115 13 2 1313 366 6910 11 9 12 7 669 13 773373711 1010 3 111010 0.0

8 8

1

5 0.1

4

0.08

Gingera Orroral Crossing Tinderry

Min. TSI

Min. TSI

0.3

1

4

8

1 4

4

4

0.02

5 5 0.00

0.2

0.4

0.6

2

12 0.04

0.8

1

1

0.06

12 13 13 1111 13 5 13 11 12 10 12 12 13 9 9 11 10 66 11 33 3 66 910 910 7 69 10 7 77 33 7 0.00

0.05

0.10

Max. TSI

0.15

0.20

Max. TSI

0.5

(c)

(d)

1

0.4

8

0.10

8 12 12 12

0.2

0.1

0.0

12 131 13 11 11 13 1 13 8 11 12 13 4 10 10 8 11 995 10 8 99555 1010 4 4 9 5 442 11 777 33 8 7733632 8 6666 222 0.0

0.2

1 1

Min. TSI

Min. TSI

0.3

8 0.05

0.00

0.4

0.6

10 10 10 9 5 25 10 5 9 9 4 2 10 4 5 9 3 77 66 59 8 3 3 2 337 2 2 6 766 7 0.00

0.8

11 8

0.05

0.10

0.15

0.20

Max. TSI

Max. TSI 0.5

1

(e) 0.4

1

0.15

11

1

1 12 0.2

0.1

0.0

0.10

11

1

Min. TSI

Min. TSI

0.3

13 13 12 13 12 11 12 11 11 4 8 11 10 44 11 13 48 77775995910 9104 3373 59 810 13 663325810 66622258 2 12

0.0

0.2

(f)

11 11

8 11 4

0.05

8

13

4 4

0.00

0.4

Max. TSI

0.6

10 5 9 10 9 9 10 8 5 5 772 9 9 10 7 6 77 5 2 8 6625322 8 10 6633 33 0.00

0.05

4 13

0.10

0.15

0.20

Max. TSI

Fig. 3. Plot of Sobol Total Sensitivity Indices (TSI) for the different target functions and Sacramento parameters for the five catchments. Minimum and maximum TSI are calculated from the four decadal values. Note (a) is for NSE, (b) zoom of (a), (c) NSElog, (d) zoom of (c), (e) F20, and (f) zoom of (e). One-to-one line is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

142

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152 0.8

Brindabella

8 8

8

8

Burbong Gingera Orroral Crossing Tinderry

0.4

0.10

3 3

0.0

2 2 2 4

0.15

Min. TSI

Min. TSI

0.6

2

(b)

8

Catchment

0.2

2

0.20

(a)

2 2 2 2 4 2 777 3 44 4 4 9 655 3 11196677 5 5 5 1 916 6 9 9 0.0

3

7 77 4 4 4 5 3

0.05

0.00

−0.05

5

6 6 5 1 6 6 7 9 119 9 7 6 1 1 9

3

4

5

5

9 −0.10

0.2

0.4

0.6

0.8

0.00

0.05

0.10

Max. TSI

0.15

0.20

Max. TSI

0.5

(c) 0.4

8

8 7

0.04

7

6

(d) 2

7

0.02

7

5 Min. TSI

0.3

Min. TSI

8 8

8 0.2

7

5

0.0

6 6 222 3 6 52 115442 6 1 551 491 96 4 5 33 994 33

−0.02

9

3 6

2 4

4

4 4 9 5 41 3 1 53 9 1 5 9

1

3

0.00

3 0.1

6

2 2 2

9

1

6

−0.04 0.0

0.2

0.4

0.00

0.6

0.05

0.10

0.15

Max. TSI

Max. TSI 0.6

(e)

0.15

7

(f)

0.5

7 8 8 8 6 6 6 8

0.3

0.2

7

9 7

0.10

7

Min. TSI

Min. TSI

0.4

4

1 2 1 1

6 9

0.1

0.0

9 4 9 8 1 24 1 14 6 9 21 21 2 3 4 3 5 2 4 3 5 33 9 555 0.0

0.2

0.00

0.4

Max. TSI

0.6

0.8

2 2 5 5 4 3 5 3 5 5 0.00

1

0.05

6

4 2 3

4 3

4 3

2

1 0.05

9

0.10

0.15

0.20

Max. TSI

Fig. 4. Plot of Sobol Total Sensitivity Indices (TSI) for the different target functions and SIMHYD parameters for the five catchments. Minimum and maximum TSI are calculated from the four decadal values. Note (a) is for NSE, (b) zoom of (a), (c) NSElog, (d) zoom of (c), (e) F20, and (f) zoom of (e). One-to-one line is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

143

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152 Table 1 TSI values using original and changed range of the IHACRES and SIMHYD parameters for the Gingera and Tinderry with the NSE objective function. IHACRES Parameter name Original 70s Gingera

f

ss sq

vs Tinderry f

ss sq

vs SIMHYD

90s

00s

70s

80s

90s

70s

80s

90s

0.38 0.05 0.55 0.50

0.80 0.03 0.30 0.29

0.06 0.05 0.62 0.71

0.10 0.04 0.75 0.62

0.11 0.05 0.71 0.58

0.31 0.05 0.62 0.54

0.47 0.04 0.36 0.47

0.74 0.03 0.30 0.28

0.69 0.04 0.35 0.33

0.92 0.02 0.20 0.18

0.29 0.05 0.51 0.57

0.49 0.03 0.53 0.44

0.41 0.05 0.57 0.48

0.81 0.03 0.30 0.28

0.23 0.04 0.45 0.65

0.43 0.04 0.49 0.50

0.35 0.05 0.55 0.53

0.76 0.03 0.33 0.32

0.23 0.01 0.44 0.48

0.45 0.01 0.34 0.42

0.45 0.01 0.39 0.43

0.96 0.01 0.11 0.16

0.04 0.01 0.53 0.55

0.11 0.01 0.49 0.58

0.15 0.01 0.52 0.58

0.90 0.01 0.12 0.17

0.42 0.01 0.37 0.39

0.68 0.01 0.23 0.27

0.74 0.01 0.27 0.28

0.99 0.01 0.13 0.17

0.23 0.01 0.46 0.45

0.47 0.01 0.34 0.36

0.48 0.01 0.38 0.38

0.95 0.01 0.11 0.15

0.21 0.01 0.44 0.50

0.43 0.01 0.32 0.41

0.45 0.01 0.37 0.42

0.94 0.01 0.13 0.21

Tinderry RISC InflCoef InfilShp SMSC IntrCoef RechCoef BaseCoef PervFrac ImpvThsh

80s

Infilshp: 0–5 90s

0.00 0.18 0.17 0.15 0.02 0.00 0.01 0.74 0.01

00s

70s

80s

90s

00s

Infilshp: 0–20

70s

80s

90s

0.00 0.18 0.26 0.13 0.07 0.03 0.08 0.71 0.00

0.00 0.21 0.17 0.08 0.08 0.01 0.12 0.76 0.00

0.00 0.21 0.18 0.12 0.12 0.02 0.11 0.80 0.00

0.00 0.19 0.14 0.12 0.14 0.03 0.11 0.82 0.00

0.00 0.00 0.15 0.14 0.15 0.10 0.17 0.15 0.02 0.15 0.00 0.04 0.01 0.04 0.77 0.76 0.01 0.01

0.00 0.16 0.03 0.11 0.01 0.00 0.01 0.78 0.01

0.00 0.18 0.04 0.15 0.02 0.01 0.02 0.81 0.01

0.00 0.01 0.15 0.09 0.03 0.03 0.15 0.12 0.03 0.23 0.01 0.07 0.02 0.06 0.83 0.92 0.01 0.01

0.00 0.01 0.01 0.20 0.21 0.18 0.37 0.38 0.34 0.06 0.09 0.10 0.04 0.05 0.06 0.02 0.02 0.03 0.10 0.12 0.11 0.67 0.68 0.70 0.00 0.00 0.00 0.00 0.17 0.13 0.13 0.01 0.00 0.01 0.73 0.01

00s

00s

70s

80s

70s

80s

ss: 10–1000

0.45 0.04 0.50 0.46

Original

RISC InflCoef InfilShp SMSC IntrCoef RechCoef BaseCoef PervFrac ImpvThsh

ss: 10–250

f: 0.5–1.7

0.25 0.05 0.48 0.61

70s Gingera

80s

f: 0.5–0.9

90s

00s

Inflcoef: 0–200 90s

00s

70s

80s

Inflcoef: 0–800 90s

00s

70s

80s

90s

0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.18 0.24 0.27 0.24 0.21 0.35 0.39 0.34 0.29 0.13 0.13 0.11 0.08 0.64 0.68 0.66 0.67 0.47 0.54 0.51 0.48 0.32 0.35 0.31 0.15 0.07 0.10 0.11 0.16 0.06 0.09 0.09 0.13 0.08 0.12 0.12 0.17 0.02 0.03 0.04 0.07 0.01 0.02 0.03 0.08 0.06 0.10 0.11 0.03 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.00 0.01 0.01 0.02 0.09 0.11 0.12 0.12 0.09 0.10 0.11 0.11 0.09 0.11 0.11 0.11 0.86 0.70 0.78 0.79 0.90 0.73 0.82 0.84 0.91 0.73 0.78 0.80 0.00 0.00 0.01 0.01 0.01 0.00 0.01 0.01 0.01 0.00 0.00 0.00

illustrative purposes. The results reflect a best case in the sense that TSI values were computed for each decade of data but only the maximum TSI value over each decade is presented here as an indication of sensitivity. The section to follow on Question 10 further explores this issue of TSI varying with respect to input forcing. The SA revealed some models to have insensitive parameters in some catchments and not in others, as reflected by the square brackets around a parameter number. Note that whether a parameter is insensitive also depends on the objective function or prediction function used. It is clear that the two more complex models can have most of their 9 and 13 parameters being insensitive. GR4J consistently has only one insensitive parameter for these three target functions while the number for IHACRES varies from one to three depending on the climate forcing of the decade and catchment. Ways of dealing with the insensitive parameters in Table 2 are discussed in the following questions. 4.3. Question 3: When can one safely fix the insensitive parameters? It is common practice to fix insensitive parameters to improve overall model identifiability, but this paradoxically assumes that those parameters are known with sufficient certainty, and that making this assumption will not impact on the prediction (Beck, 1987). Insensitivity indicates that the parameters have little impact on objective function values, or on the predictions that are made using the model. The first case indicates how well each parameter can be calibrated using a particular objective function, while the second indicates how important the parameter is in terms of the purpose of the modeling exercise. Ideally, parameters that are important need to be well determined, while parameters that are not important can have considerable uncertainty without impacting on the conclusions of the study. If there are any parameters that are important in terms of the conclusions being made,

0.01 0.16 0.36 0.11 0.01 0.00 0.01 0.60 0.02

0.00 0.19 0.43 0.14 0.01 0.00 0.01 0.69 0.01

0.00 0.01 0.16 0.12 0.41 0.64 0.15 0.09 0.01 0.14 0.00 0.04 0.01 0.06 0.72 1.00 0.01 0.01

0.00 0.24 0.18 0.10 0.01 0.00 0.01 0.66 0.02

0.00 0.29 0.25 0.14 0.02 0.00 0.01 0.72 0.01

00s

0.00 0.00 0.25 0.25 0.23 0.36 0.15 0.12 0.02 0.19 0.01 0.06 0.01 0.06 0.75 0.94 0.01 0.01

0.00 0.12 0.08 0.15 0.01 0.00 0.01 0.77 0.01

0.00 0.13 0.11 0.16 0.02 0.01 0.02 0.79 0.01

00s 0.00 0.14 0.23 0.17 0.15 0.04 0.09 0.85 0.00

0.00 0.01 0.11 0.08 0.09 0.10 0.17 0.15 0.02 0.22 0.01 0.06 0.02 0.06 0.81 0.91 0.01 0.01

but do not contribute significantly to the objective function, then consideration needs to be given to modifying the objective function, adopting a completely different objective function, or employing multiple objective functions. If the conclusions are not sufficiently dependent on the value of a parameter (uncertainty in the conclusion is considerably larger than the uncertainty due to the particular parameter), then from the point of view of the conclusions being made, the parameter value can be set to a ‘‘reasonable’’ value that is based on other information. In Fig. 1, the ss (no. 2) parameter of the IHACRES model has maximum TSI values for the objective functions NSE and NSElog lower than the threshold 0.2, so it may be difficult to identify. The maximum TSI value for the F20 is also lower than 0.2 (except for the Burbong catchment), so ss has a small effect on the frequency of low flows. As both conditions hold, that is a parameter is insensitive for both the objective and prediction functions, ss can safely be fixed for calibration (excepting Burbong). Similarly, in Fig. 2 with the GR4J model, the X4 (no. 4) parameter in Gingera is insensitive with respect to the NSE objective function and for the F20 prediction function in all decades, and could therefore be safely fixed to improve identifiability. In Figs. 3 and 4, the five parameters of the Sacramento model (UZK (no. 3), ZPERC (no. 6), REXP (no. 7), LZFSM (no. 9) and LZFPM (no. 10)) and one parameter of the SIMHYD model (RISC (no. 1)) are always insensitive in relation to all our target functions so they can be safely fixed if those functions capture the model purpose. Note that the parameter SMSC (no. 4) of the SIMHYD model is excluded because it is located near the threshold 0.2 so care needs to be taken. However the LZPK (no. 12) parameter of Sacramento and the BaseCoef (no. 7) parameter of SIMHYD are insensitive with respect to the NSE objective function but they are sensitive for the F20 prediction function for five of the catchments. Consequently, it may not be desirable to fix these parameters of our two more

144

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

complex models when calibrating with NSE for these catchments; while they are insensitive with respect to the objective function statistics they could lead to serious prediction errors in the prediction function to be reported to model users. Other measures to improve identifiability need to be considered in this case, such as invoking another objective function. 4.4. Question 4: Which objective function(s) will be most sensitive to help identify a particular parameter? The sensitivity of parameters depends on the objective functions used, hence different objective functions could help in identifying parameters. van Werkhoven et al. (2008, 2009) found that different objective functions affect the sensitivity of different parameters of the Sacramento model. They used four different objective functions: root mean squared error (RMSE) to emphasize the high flow fitting; transformed root mean squared error (TRMSE) to emphasize the low flow fitting; runoff coefficient error (ROCE) to capture the water balance; and slope of the flow duration curve error (SFDCE) to capture the distribution of mid-level flow. For the IHACRES model (see Fig. 1), the ss (no. 2) parameter, that represents the time constant for the slow flow store, is insensitive with respect to the NSE objective function, largely because the objective function is focused on fitting high flows. But it becomes more sensitive with respect to the NSElog objective function because that objective function is more focused on low flow properties. For the GR4J model (see Fig. 2), the X4 (no. 4) parameter, that represents the time base of the unit hydrograph, is insensitive with respect to the NSElog objective function. But it becomes more sensitive with respect to the NSE objective function because the time base is more sensitive towards capturing the high flow peak. Not surprisingly, therefore, the parameters in a model that are related to low flows are insensitive to a high flow objective function but become sensitive with respect to low flow objective function and vice versa. Similar results apply with the two more complex models. The Sacramento model (see Fig. 3) has insensitive parameters LZPK (no. 12) and PFREE (no. 13)) with the NSE objective function but they become more sensitive with the NSElog objective function. For the SIMHYD model (see Fig. 4), BaseCoef (no. 7) is insensitive with the NSE objective function but it become more sensitive with the NSElog objective function. Therefore sensitivity analysis can give valuable insight into the relationship between the impact of a parameter on the model output and the behavior captured by an objective function. For example, an objective function that measures model performance at high flows can be used to calibrate parameters controlling high flows. Sensitivity analysis can be used to design objective functions that target particular model parameters leading to potentially a multi-criteria calibration of the model, depending on the potential of sensitivity analysis to create subsets of parameters that influence different aspects of the catchment response. 4.5. Question 5: If a hybrid objective function is used, are there fewer insensitive model parameters? From Fig. 1 with the IHACRES model, it can be seen that the hybrid objective function, 0.5NSE + 0.5NSElog, increases the sensi-

tivity of the f (no. 1) parameter compared to the NSE objective function alone, due to equal weight being given to the NSElog objective function. The f (no. 1) parameter controls the amount of actual evaporation and is related to the drying of the catchment, so it affects the response of the catchment in drier conditions; the NSElog objective function reflects the fitting of low flows more. However, the sensitivities of the sq (no. 3) and vs (no. 4) parameters that are more related to high flow responses are reduced with use of the hybrid objective function compared to the NSE objective function. (Note that the vs (no. 4) parameter, the fractional volume for slow flow, is related to high flow because the vs is calculated as: vs = 1  vq). In fact, the Tinderry catchment has two more insensitive parameters, and the Burbong catchment has one more insensitive parameter, with use of the hybrid objective function compared to the NSE objective function. It is evident that a hybrid objective function can increase the sensitivity of some parameters while decreasing the sensitivity of other parameters. Therefore, a hybrid objective function is not a panacea for converting all parameters to be more sensitive, and a complementary strategy would be for modelers to select the insensitive parameters that need to be fixed, taking modeling purpose into account. 4.6. Question 6: Does it matter whether one selects Morris or Sobol? Do they give similar results? Using the easy and simple implementation of the code in Appendix C, Fig. 5 shows the overall sensitivity of each parameter in the case of: the Morris method, the NSE objective function, the Gingera catchment in the 1980s, and all four rainfall–runoff models. The numbers in the figure identify each parameter, as shown in Table A1. The result of applying the Morris method is similar to that of Sobol (see TSI values in Figs. 1–4) but there is a slight difference in the ordering of the parameter sensitivities. For example, the order of sensitivity for the case of Sacramento by the Morris method is parameter no. 2 > no. 8 > no. 5 while Sobol TSI values yield a slightly different order: no. 8 > no. 5 > no. 2. As Campolongo and Saltelli (1997) have stated there can be a discrepancy in the identification of parameters between the Morris and Sobol algorithms. They reported that several parameters that are well identified by Morris may not be well identified by Sobol and vice versa. Both Sobol and Morris methods sometimes produce ‘‘type I error’’ (incorrectly classifying a sensitive parameter as being insensitive), but they do not produce ‘‘type II error’’ (incorrectly classifying an insensitive parameter as being sensitive). For this paper, we preferred Sobol because, as described in its description, it provides ratio-scaled quantified values and computational time was not a constraining issue. These results however support the observation that Morris may be useful for screening by reducing computational time, yet provides a similar, though subjectively interpreted, identification of insensitive parameters. 4.7. Question 7: To obtain reliable estimates of sensitivity, how long does the data period of analysis need to be? Figs. 6 and 7 show Box plots of Sobol TSI values for different durations of years from one to 40 years to discover what constitutes a reasonable duration for sensitivity analysis. A whole

Table 2 Insensitive parameters (max TSI < 0.2) with respect to NSE, NSElog and F20 target functions.a

a

Model

NSE

NSElog

F20

IHACRES GR4J Sacramento SIMHYD

ss (no. 2 in Table A1)

ss (no. 2), [sq (no. 3)], [vs (no. 4)]

[X4 (no. 4)] All except for nos. [1], [2], [4], [5], 8 All except for nos. [2], [3], [5], 8

X4 (no. 4) All except for nos. 1, [4], [11], 12, 13 All except for nos. [6], 7, 8, [9]

[ss (no. 2)], sq (no. 3) X4 (no. 4) All except for nos. 1, [4], 12, [13] All except for nos. [6], 7, 8, [9]

The square bracket [ ] around a parameter indicates that it is insensitive for some catchments but not others.

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

40-year period is divided into shorter periods of one, two, five, 10, or 40 years, hence the results for two (five) years of duration correspond to SA being applied on 20 (eight) consecutive periods (and so on). This corresponds to the process followed by Kim et al. (2011). The repeated periods address climate variation in the historical record, but the single catchment used obviously limits the breadth of conditions considered, and therefore limit the transferability of these results. The figures are therefore illustrative results for the IHACRES and SIMHYD models applied to the Gingera catchment, which displays a broad range of climatic conditions (see Table A2). The TSI values corresponding to our maximum possible 40-year duration is plotted using a dot-dash gray line across the other duration Box plots. From the figures, one and two year durations have quite a large range of TSI variation for most of the parameters, and in these cases one must analyse sensitivity for a five to 10 year period in order to achieve a stable and similar range of TSI variation. Indeed five years may be the minimum duration of daily rainfall data needed to obtain stable TSI values from sensitivity analysis with IHACRES and SIMHYD models, at least with the NSE objective function for the Gingera catchment. Note that the results could be affected by changing the parameter ranges investigated, as seen for the results from Question 1. 4.8. Question 8: Which rainfall–runoff models tend to be more difficult to identify? From Fig. 1, the ss (no. 2) parameter of the IHACRES model is insensitive with respect to the NSE objective function but the TSI value of this parameter increases up to 0.2 when the NSElog objective function is invoked. The X4 (no. 4) parameter of the GR4J model (Fig. 2) has low TSI value with the NSElog objective function but has a TSI value over the threshold with the NSE objective function. The presence of insensitive parameters with a certain objective function does not imply that poor identification cannot be avoided. Those parameters may be sensitive with respect to other objective functions. However if the parameters are insensitive with respect to all the objective functions of interest then these parameters could have identification problems. In Fig. 3, Sacramento has 8 out of 13 parameters insensitive for both NSE (nos. 3, 6, 7, 9, 10, 11, 12, 13) and NSElog (nos. 2, 3, 5, 6, 7, 8, 9, 10). In Fig. 4, SIMHYD has 5 out of 9 parameters insensitive for both NSE (nos. 1, 4, 6, 7, 9) and NSElog (nos. 1, 2, 3, 4, 5). Of these, five parameters of the Sacramento model, UZK (no. 3), ZPERC (no. 6), REXP (no. 7), LZFSM (no. 9) and LZFPM (no. 10) are always insensitive with respect to both objective functions as we discussed in Question 3 (see the expanded plot in Fig. 3 where parameter nos. 3, 6 and 7 are clustered in the bottom left corner, and nos. 9 and 10 also have maximum and minimum TSI lower than 0.1). Parameter RISC (no. 1) of the SIMHYD model is insensitive with both objective functions (note that the parameter SMSC (no. 4) is excluded because it is located near the threshold 0.2 as we mentioned early in Question 3). From this point of view, the Sacramento and SIMHYD models could be regarded as overparameterised models which may be difficult to identify. They have a large number of insensitive parameters with each objective function, and some of these parameters are insensitive with all our objective functions. Additionally we could examine model structure issues by sensitivity analysis. From Fig. 1, the sensitivity of the ss (no. 2) parameter in IHACRES for the drier catchments (Tinderry and Burbong) is lower than the threshold of 0.2 for both maximum and minimum TSI with respect to the NSElog objective function. But this parameter, representing the time constant for the slow flow store, is not sufficiently sensitive with respect to the NSElog (low flow) objective function. This means that the structure of IHACRES may be difficult to calibrate when modeling low flows in dry catchments. For

145

improvement of low flow fitting, additional parameters for groundwater modeling, for example, could be one of the solutions. However for the case of the GR4J model (Fig. 2) with the Burbong catchment and the NSElog objective function, all parameters except for X4 (no. 4) have TSI values above the threshold. The parameter X4 (no. 4) is more sensitive for capturing the high flow peak as explained in Question 4. Therefore the parameters of the GR4J model work well for low flow modeling. For the Tinderry catchment, the parameter X1 (no. 1) has TSI values slightly below the threshold in terms of minimum TSI but this is due to there not being enough information in the data for the very dry period of the 2000s. Hence the GR4J model seems more suitable than the IHACRES model for modeling low flows in dry catchments, and the IHACRES model may benefit from additional parameters related to groundwater interactions with streamflow. These results therefore support prior observations that overparameterised model structures are problematic, and simpler conceptual model structures may be preferable (Jakeman and Hornberger, 1993). At the same time the best model structure for a given catchment or dataset may vary depending on its characteristics (Perrin et al., 2001). 4.9. Question 9: Does the sensitivity of parameters vary between catchments? It was found by van Griensven et al. (2006) that the sensitivity of parameters of a SWAT model varies between catchments. For a more informed parameter sensitivity analysis and identification, examination of individual catchment datasets is needed. Let us examine this finding for the case of the IHACRES rainfall–runoff model (see Fig. 1). For the two wettest catchments (Brindabella and Gingera) with the NSE (high flow) objective function, the sq (no. 3) and vs (no. 4) parameters (that are related to high flow) have TSI values above the threshold 0.2 in terms of the both maximum and minimum TSI. This implies that the data of the wetter catchments have enough information over the four different decadal datasets to identify the high flow parameters of the IHACRES model, given their effect on the high flow objective function. With respect to the data from the dry catchments (Tinderry and Burbong) and the NSE objective function, the sq (no. 3) and vs (no. 4) parameters have their maximum TSI value above the threshold but have their minimum TSI value under the threshold. This implies that there is insufficient high flow information in some decadal datasets in dry catchments. Our results therefore support the observation of van Griensven et al. (2006) that sensitivities of parameters in a rainfall–runoff model structure are specific to a site, and cannot be assumed from previous work in other catchments. 4.10. Question 10: If a different data period is used, can there be fewer insensitive model parameters? The variability of parameter sensitivity of the four different rainfall–runoff models with four different decadal datasets is examined here using the maximum and minimum TSI values from the SA undertaken (Figs. 1–4). Each point is assigned coordinates using maximum TSI (horizontal axis) and minimum TSI (vertical axis). If the maximum and minimum TSI values are equal for a parameter, that is they are on or near the diagonal line, that parameter has the same TSI for all four decades. If some parameters of a model are insensitive for all data periods then the parameters cannot be identified by the calibration process. However if the parameters are sensitive for a particular data period then these parameters could be identifiable with that data, assuming that the same parameter values are appropriate across periods. Note that periods where different hydrological

146

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

3

3

(b)

1.2

1.0

(a)

0.6

0.2

0.4

0.4

1

2

0.8

σ

4

0.6

0.8

1.0

1

2 0.3

0.4

0.5

0.6

1.4

0.2

4 0.2

2

0.4

0.5 8

(d)

1.0

0.8

1.2

(c)

0.3

σ

0.8

0.6

2

0.4

0.6

4 8

0.4

5 12

9 10 11 13 3 7 6

0.0

0.1

5 1 4

0.2

7

0.2

0.2 0.0

3

6 1 9 0.3 μ

0.4

0.5

*

0.0

0.1

0.2

0.3

0.4 μ

0.5

0.6

*

Fig. 5. Morris sensitivity index for four different rainfall–runoff models using the 1980s decadal period, and for the Gingera catchment with the NSE objective function. Note (a) IHACRES, (b) GR4J, (c) Sacramento, and (d) SIMHYD. Numbers correspond to the parameters described in Table A1 in Appendix A.

f

1.0

0.25

0.8

0.20

0.6

0.15

0.4

tau_s

0.10

0.2

yr

0.05

1

TSI

0.0 tau_q

2

v_s

5

0.6

10

0.6

0.5 0.4

0.4

0.3 0.2

0.2 0.1 1

2

5

10

1

Duration (Year)

2

5

10

Fig. 6. Variation of TSI with different data period durations (in years) for the IHACRES rainfall–runoff model with Gingera catchment and the NSE objective function.

behavior is observed may not have the same parameter values, and calibrating on those periods would not be appropriate. For example, the climate forcing and hence response for the 2000s in the Tinderry catchment differ significantly from that in the other three decades. There will however be occasions where variation in sensitivity between periods might be used advantageously to improve identifiability. Note that this differs from Question 7, where the variation of TSI value with duration of data is undesirable because

the instability of TSI value results from there not being enough information in the data due to improper selection of period length. With a duration of 10 years pre-selected (as described in Section 3, and greater than the minimum of five years for sensitivity analysis identified in Question 7), we now evaluate the variation of TSI value across the four decades, including wet and dry periods. Parameters that are insensitive in one decade might become sensitive in another. This result is obviously highly dependent on the catchment’s temporal characteristics and data available. Using decadal data, there is no such opportunity for changing the calibration data period with the SIMHYD model calibrated with the NSE objective function for the Tinderry catchment (see Fig. 4a). All parameters except for no. 8 are insensitive in terms of both minimum and maximum TSI values, so these parameters are insensitive regardless of the data period used. The nos. 2 and 4 parameters especially are near the diagonal line that represents the maximum and minimum TSI values, indicating similarity for all four decades. This model would not have fewer insensitive parameters if a different data period is used. For the parameters of the IHACRES model, again with the NSE objective function and the Tinderry catchment (see Fig. 1a), sq (no. 3) and vs (no. 4) parameters are insensitive in terms of minimum TSI. For the parameters of the GR4J model with the NSE objective function and Tinderry catchment (see Fig. 2a), the X4 (no. 4) parameter is insensitive in terms of minimum TSI. In both cases, this occurs only when using 2000s data (we do not show the TSI value for all decadal periods), which is known to have different hydrological behavior. Calibrating on other periods to improve identifiability could therefore result in poor performance in the 2000s. But the parameters are sensitive in all other periods, so there is no need to improve identifiability for those other decadal periods.

147

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

RISC

InflCoef

0.01

0.25

0.00

0.20

−0.01

InfilShp 0.6 0.4

0.15

−0.02

0.2

0.10

−0.03

0.05

0.0

−0.04 SMSC

0.4

0.20

RechCoef yr

0.06

0.3

0.15

TSI

IntrCoef

1

0.04

0.10

0.2

0.05

0.1

2

0.02

5

0.00

10

0.0 BaseCoef

PervFrac

0.20

ImpvThsh

0.8

0.00

0.7

−0.02

0.6

−0.04

0.15 0.10 0.05 1

2

5

10

1

2

5

10

1

2

5

10

Duration (Year) Fig. 7. Variation of TSI with different data period durations (in years) for the SIMHYD rainfall–runoff model with Gingera catchment and the NSE objective function.

Table 3 Answers to the 10 questions. Questions

Answers

Question 1: Is the result of a sensitivity analysis affected by the parameter ranges selected?

Yes, changing the range of sensitive parameters can make the parameters more sensitive or insensitive and should therefore include only plausible parameter sets, as determined by parameter uncertainty. Using multiple parameter ranges can also help characterize uncertainty in the sensitivities Each target function is insensitive to some, often different, parameters, particularly for those models with more than a handful of parameters Insensitive parameters can only be fixed if the prediction function of interest is also insensitive to that parameter, which is the case with several of these models The objective function that relates to the same purpose or aspect that the parameter captures

Question 2: Is the selected target function(s) insensitive to changes in any parameters? Question 3: When can one safely fix the insensitive parameters? Question 4: Which objective function(s) will be most sensitive to help identify a particular parameter? Question 5: If a hybrid objective function is used, are there fewer insensitive model parameters? Question 6: Does it matter whether one selects Morris or Sobol? Do they give similar results? Question 7: To obtain reliable estimates of sensitivity, how long does the data period of analysis need to be? Question 8: Which rainfall–runoff models tend to be more difficult to identify?

Question 9: Does the sensitivity of parameters vary between catchments? Question 10: If a different data period is used, can there be fewer insensitive model parameters?

No, some parameters are sensitive with a hybrid objective function but other parameters become less sensitive due to the reduced weight of the contributing objective functions Both of them give similar results except for a slightly different ordering of parameter sensitivity, but the discrepancy of order does not affect the identification of sensitive parameters Five years is the minimum period required for securing sensitivity analysis results that are stable This SA supports prior observations that model structures are problematic if they have many insensitive parameters with any objective functions, e.g. Sacramento and SIMHYD, and simpler conceptual model structures may be preferable. At the same time the best model structure for a given catchment or dataset may vary depending on its characteristics, e.g. GR4J for low flows in dry catchments rather than IHACRES Yes, and so results from previous uses of a model structure cannot be assumed in other case studies Yes. Depending on the model structure and catchment, this SA supports the suggested practice of selecting different data periods to overcome identifiability problems

For the Sacramento model, examining the sensitivity of NSE in Tinderry reveals that some parameters are insensitive in all periods (maximum TSI value is lower than 0.2), as with the SIMHYD model, and some are only insensitive in the 2000s, as with the IHACRES and GR4J models. Parameter no. 1 (UZTWM) is insensitive in the 1970s and becomes sensitive in the 1980s and 1990s. Calibrating in those periods would therefore improve the identifiability of this parameter, though other parameters remain insensitive, and possibly difficult to identify. In this context, Wagener et al. (2003) developed an approach of dynamic identifiability analysis to infer, for a

given objective function, periods in a dataset that provided high information content to identify parameters. 5. Discussion and conclusions This paper has presented results of sensitivity analysis using the Hydromad and Sensitivity packages in R. A sensitivity analysis can be performed with just a few lines of code after defining a rainfall– runoff model structure. We provided some answers to 10 possible questions for the sensitivity analysis using four different target

148

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

Table A1 Description of parameters of the four models. Parameter name

Parameter number

Unit

Range

Description

IHACRES-CMD f e d ss (tau_s) sq (tau_q) vs (v_s)

1 – – 2 3 4

[–] [–] [mm] [day] [day] [–]

0.5–1.3 1 (fixed) 200 (fixed) 10–500 0–10 0–1

CMD stress threshold as a proportion of d Temperature to potential evapotranspiration (PET) conversion factor CMD threshold for producing flow Time constant for slow flow store Time constant for quick flow store Fractional volume for slow flow

GR4J X1 X2 X3 X4

1 2 3 4

[mm] [mm] [mm] [day]

50–5000 15 to 4 10–1300 0.5–5

Maximum capacity of the production store Groundwater exchange coefficient One day ahead maximum capacity of the routing store Time base of unit hydrograph UH1

Sacramento UZTWM UZFWM UZK PCTIM ADIMP ZPERC REXP LZTWM LZFSM LZFPM LZSK LZPK PFREE SIDE RSERV RIVA

1 2 3 4 5 6 7 8 9 10 11 12 13 – – –

[mm] [mm] [1/day] [–] [–] [–] [–] [mm] [mm] [mm] [1/day] [1/day] [–] [–] [–] [–]

1–150 1–150 0.1–0.5 0.000001–0.1 0–0.4 1–250 0–5 1–500 1–1000 1–1000 0.01–0.25 0.0001–0.25 0–0.6 0.0 (fixed) 0.3 (fixed) 0.0 (fixed)

Upper zone tension water maximum capacity Upper zone free water maximum capacity Upper zone free water lateral depletion rate Fraction of the impervious area Fraction of the additional impervious area Maximum percolation rate coefficient Exponent of the percolation equation Lower zone tension water maximum capacity Lower zone supplementary free water maximum capacity Lower zone primary free water maximum capacity Lower zone supplementary free water depletion rate Lower zone primary free water depletion rate Direct percolation fraction from upper to lower zone free water storage Fraction of base flow that is draining to areas other than the observed channel Fraction of the lower zone free water that is unavailable for transpiration purposes Fraction of the riparian vegetation area

SIMHYD RISC InflCoef InfilShp SMSC IntrCoef RechCoef BaseCoef PervFrac ImpvThsh

1 2 3 4 5 6 7 8 9

[mm] [–] [–] [mm] [–] [–] [–] [–] [mm]

0–5 0–400 0–10 1–500 0–1 0–1 0–1 0–1 0–5

Rainfall interception store capacity Infiltration coefficient Infiltration shape Soil moisture store capacity Interflow coefficient Recharge coefficient Baseflow coefficient Pervious fraction Impervious threshold

Table A2 Characteristics of the five study catchments.a

a b c

ID

Area (km2)

Mean annual P (mm)

Mean daily Q (mm)

Mean annual ROC (Q/P)

Mean decadal BI (BQ/TQ)

Mean decadal Pb (mm)

Brindabella Gingera Orroral Crossing Tinderry Burbong 1970sc 1980s 1990s 2000s

427 148 90 490 505 – – – –

1127 985 885 716 664 954 913 901 734

0.96 0.77 0.33 0.35 0.19 0.72 0.56 0.53 0.26

0.31 0.29 0.13 0.18 0.11 0.27 0.21 0.20 0.11

0.54 0.63 0.63 0.21 0.19 – – – –

1175/1179/1193/962 1061/1013/1030/837 975/929/884/749 841/759/693/572 716/684/706/549 – – – –

P: precipitation; Q: streamflow; ROC: runoff coefficient; BI: baseflow index calculated by IHACRES. The four values of mean decadal P are for the 1970s/1980s/1990s/2000s. The decadal period values are averaged values over the five catchments.

functions, five catchments and four rainfall–runoff models. Some answers to the 10 questions are summarized in Table 3. While this paper focused on rainfall–runoff modeling, the method and many of the questions are generalisable to other models. While the code in Appendix C is specific to the Hydromad rainfall–runoff modeling package, the ‘Sensitivity’ R package (Pujol et al., 2012) can equally easily be used with other models written in R. The issue in Question 1 is likely to be experienced by other models having parameters with nonlinear effects; Question 2–5 might be useful for any model requiring parameter estimation from observations; Question 6 applies to any model; the issues of over-parameterisation and model-dependent results in Ques-

tions 8 and 9 may apply in many other models; Question 7 and 10 apply to time series, but perhaps also to determining sample size more generally. Viewed at a conceptual level, the questions we attempted to answer highlight two key issues in modeling: identifiability and uncertainty in SA. Identifiability of models requiring parameter estimation from observations can be improved by, amongst other solutions, changing objective functions, fixing parameters, changing model structures or changing data periods. This paper illustrates how sensitivity analysis can play a key role in understanding the potential effectiveness of each of these measures. Uncertainty is demonstrated by results in SA changing due to the

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

selection of parameter ranges, data sample size and SA method. When using SA, a modeler should be aware of these potential limitations, and interpret results accordingly. This paper therefore brings out the simplicities and complexities of applying sensitivity analysis. Note that parameter (non)stationarity to climate change is beyond the scope of this paper. We are pointing out that this is another issue which needs to be addressed. It should be noted that this cannot be done from the view point of SA and it needs careful consideration of the impacts of climate change on the catchment (e.g. stocking densities, vegetation distributions, and evapotranspiration rates, etc.). Performing an SA is technically straightforward in R with the Sensitivity and Hydromad packages. The four questions from Question 2 to 5 demonstrate the simplicity of using SA to inform whether parameters are insensitive, which parameters can be fixed, and whether other objective functions may improve sensitivity. The remainder of the paper however demonstrates how strongly sensitivities are dependent on parameter ranges, length of data period, catchments, model structures and climatic conditions. While a modeler can easily use SA to help assess and improve their model in a particular case study, great care must be taken in drawing generalisations or trying to apply generalisations to new situations. Acknowledgements The first author was funded through a scholarship from CSIRO Land and Water, the Australian National University, and National Centre for Groundwater Research and Training (NCGRT). The second, third and fourth authors were supported by NCGRT funding within the uncertainty research program. We gratefully acknowledge Francis Chiew for supplying the source code of SIMHYD model, Dario Mavec for technical support in generating the figures, and Nathalie Saint-Geours for valuable discussion. We also gratefully acknowledge ACTEW Corporation for supplying the rainfall data covering the five study catchments, and streamflow data for the Gingera and Tinderry catchments. We also gratefully acknowledge the ACT government for supplying the streamflow data for the Burbong and Orroral Crossing catchments, and the NSW government for supplying the streamflow data for the Brindabella catchment. Appendix A. The hydrological models Conceptual models are by far the most common model type used in hydrological modeling applications. Four well-known and documented conceptual models were adopted for the investigations. Conceptual models represent flows through idealized catchment stores that are configured in parallel and/or series. Store parameters are constant in space and time for the catchment being modeled. The four models were selected because of their varying complexity, from four to 13 parameters. The hydrological models and the evolutionary algorithms for calibration are available in Hydrological Model Assessment and Development (Hydromad), an open source, R-based, software package (Andrews et al., 2011), except for the SIMHYD model which is called from Hydromad. This study uses Hydromad for all of the modeling undertaken. A description of the parameters in each model, including their range, is given in Table A1. All of the models use rainfall and potential evapotranspiration time series as input data. The latter series is calculated from pan evaporation data. A.1. IHACRES-CMD with four parameters The IHACRES rainfall–runoff model runs on any time step, having been used on six-minute time steps in very small catchments (Hansen et al., 1997) to monthly in the Thames catchment (Little-

149

wood and Marsh, 1996). Since the model’s development (Jakeman et al., 1990; Jakeman and Hornberger, 1993) in the early 1990s it has been adopted in various studies (e.g. Viney et al., 2009; Vaze et al., 2010; Petheram et al., 2012; Post et al., 2012; Silberstein et al., 2012). There are several versions available, their use depending on the type of catchment and application (Ye et al., 1997; Evans and Jakeman, 1998; Croke and Jakeman, 2004). This study uses the catchment moisture deficit (CMD) version of IHACRES (Croke and Jakeman, 2004) where rainfall is partitioned into drainage and evapotranspiration, and changes in catchment moisture store are accounted for at each time step. Effective rainfall is generated from a nonlinear function of the raw rainfall and the CMD output from the accounting equation. Separate linear unit hydrographs in a parallel storage configuration convert effective rainfall into quick and slow streamflow which can be summed at each time step for total streamflow. Other configurations of storages in parallel and/or series are possible but two parallel storages were identified as preferable for the five catchments here. Such a version can have up to six parameters notionally but the parameter e is fixed here as unity because potential evapotranspiration data is used, not temperature, and the parameter d is fixed here as 200 based on a study of Croke and Jakeman (2004) (see Table A1). The range of the parameter f (see Table A1) was selected based on the study of Croke and Jakeman (2004). The parameter vs (see Table A1) represents the ratio of the volume of slow flow to total flow, so it has a range of zero to one. The range of the parameters ss and sq (see Table A1) were set after considering the size of the catchments used in this study and therefore their likely response times. A.2. GR4J with four parameters The GR4J model (Perrin, 2000), a daily lumped rainfall–runoff model with four parameters, is a modified version of the GR3J model that was developed by Edijatno et al. (1999). This model uses two stores and two unit hydrographs. The production store, for storage in the surface of soil, has processes for storing rainfall, evapotranspiration and percolation. The routing store has one effective rainfall as input that is routed by a unit hydrograph. One unit hydrograph routes 90% of effective rainfall (slow flow that infiltrates into the ground) and the other unit hydrograph routes the remaining 10% of effective rainfall (quick flow that flows on the surface of the soil) (Le Moine et al., 2008). More details of the GR4J model are provided in Perrin et al. (2003). Various studies using this model include Le Moine et al. (2008), Thyer et al. (2009), Moussu et al. (2011), Andréassian et al. (2012) and Pushpalatha et al. (2012). The range of parameters of GR4J model in this study (Table A1) is different to the one in the paper of Perrin et al. (2003) because the latter could be too restricted for a number of catchments. Therefore larger intervals should be used to avoid optimized values being artificially stuck at the boundary of the domain, which may limit modeling efficiency. The enlarged ranges of parameters were suggested by Perrin (personal communication, 2011) and correspond to the 98% confidence intervals on a set of about 1000 French catchments. A.3. Sacramento with 13 parameters The Sacramento Soil Moisture Accounting (SAC-SMA) model was developed by Burnash et al. (1973). It is currently used by the National Weather Service River Forecast System (NWSRFS) Centers in the United States to perform not only real-time river and flood forecasts but also long term predictions. The version used here has 13 parameters whose ranges are taken from Blasone et al. (2008) (see Table A1). The standard version of the Sacramento model has 16 parameters and the three parameters SIDE, RSERV

150

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

and RIVA are fixed in this study based on the study of Peck (1976). As with the other three models, this model has also been used extensively in research studies (e.g. van Werkhoven et al., 2009; Vaze et al., 2010; Petheram et al., 2012; Post et al., 2012; Silberstein et al., 2012). The Sacramento model uses soil moisture accounting for water balance simulation by distributing soil moisture in various depths of the soil and is represented as a network of interconnected tanks. The soil moisture store is filled by rainfall, and reduced by evaporation and flow out from the store. The excess rainfall becomes runoff that is transformed by a unit hydrograph. The lateral flow from the soil moisture store is superimposed on this runoff and turns into streamflow (Podger, 2004). The Sacramento model is constituted by an upper and a lower zone, each including tension and free water storages that interact with each other. The algorithm of this model consists of five different runoff components: (1) direct runoff from permanent and temporary impervious areas, (2) surface runoff of excess rainfall and interflow when upper zone stores are full, (3) interflow of free water, (4) supplemental base flow, and (5) primary base flow. The free water defines the water that is free to move through the soil. The tension water represents the water that is held in position by the molecular attraction between soil particles and water molecules (Burnash, 1995). A.4. SIMHYD with nine parameters The SIMHYD model (Chiew et al., 2002) is a lumped conceptual daily rainfall–runoff model. This model is a simplified version of the HYDROLOG model that was developed in the early 1970s (Porter, 1972; Porter and McMahon, 1975) and it has been used in many studies (e.g. Zhang and Chiew, 2009; Chiew et al., 2010; Petheram et al., 2012; Post et al., 2012; Silberstein et al., 2012). This study uses a version of the SIMHYD model with nine parameters whose ranges are taken from Podger (2004). For a brief description of the SIMHYD model see Podger (2004). Daily rainfall charges the interception store which is emptied each day by evapotranspiration. The excess rainfall is used by an infiltration function that determines the infiltration capacity. The excess rainfall that exceeds the infiltration capacity turns into infiltration excess runoff. The infiltrated water is used by a soil moisture function that diverts the water to the three components: the stream (interflow), the groundwater store (recharge) and the soil moisture store. The interflow is firstly estimated by a linear function of the soil wetness (soil moisture level divided by soil moisture capacity). Then groundwater recharge is estimated also by a linear function of the soil wetness. The remaining water flows into the soil moisture store. The evapotanspiration from the soil moisture store is estimated by a linear function of the soil wetness, and it cannot exceed potential evapotranspiration. The soil moisture store has a finite capacity, hence the overflow transfers into the groundwater store. Baseflow from the groundwater store is calculated using a linear recession equation from the store. The total

runoff is the sum of impervious runoff, infiltration excess runoff, interflow (and saturation excess runoff) and baseflow. Appendix B. Study area The five catchments selected for the investigations are the Goodradigbee River at Brindabella (here Brindabella), Cotter River at Gingera (Gingera), Orroral River at Crossing (Orroral Crossing), Queanbeyan River at Tinderry (Tinderry) and Molonglo River at Burbong (Burbong) in the Australian Capital Territory (ACT). The range of watershed sizes is from 90 to 505 km2 (see Table A2): Gingera and Orroral Crossing are around 100 km2 and the other three catchments are larger than 400 km2. All the catchments are in the mountainous region except for Burbong and the western three catchments are at a higher elevation than the Tinderry catchment. There is a marked difference in rainfall amount from the western two catchments, Brindabella and Gingera with a higher rainfall regime (985 – 1127 mm/yr of mean annual rainfall), the intermediate catchment Orroral Crossing (885 mm/yr of mean annual rainfall), to the eastern two catchments, Tinderry and Burbong with a lower rainfall regime (664–716 mm/yr of mean annual rainfall). The runoff responses to rainfall are markedly different also. For example the mean annual runoff coefficient (proportion of rainfall) and mean decadal baseflow indices (ratio of baseflow to total flow in the streamflow) for the five catchments are: Brindabella (0.31, 0.54), Gingera (0.29, 0.63), Orroral Crossing (0.13, 0.63), Tinderry (0.18, 0.21) and Burbong (0.11, 0.19). The decadal baseflow indices are from calibrated flow using the IHACRES model. The differences in rainfall amount between the 1970s and the 2000s are significant. The mean annual rainfall for the 1970s and the 2000s are 954 mm, 913 mm, 901 mm and 734 mm respectively. Hence the amount of rainfall in the 2000s represents a 23% reduction compared to the 1970s. Appendix C. The use of SA algorithms using the R Sensitivity package R is a free software environment for statistical computing and graphics. The Hydromad and Sensitivity packages make it easy to perform sensitivity analyses of conceptual rainfall–runoff models by working within the flexibility of a full-featured programming environment for power-users while hiding the complexity of the underlying model and algorithms for new users. R is available from http://www.r-project.org/. Instructions for installing Hydromad can be found at http://hydromad.catchment.org/. This appendix presents the R code used to perform a sensitivity analysis. install.packages(‘‘sensitivity’’)

library(sensitivity) library(hydromad)

 Install Sensitivity package from the Comprehensive R Archive Network (CRAN)  Load Sensitivity package  Load Hydromad package

Define a conceptual rainfall–runoff model structure to use with an existing dataset in ‘zoo’ format: model.str <- hydromad(obs, sma = ‘‘cmd’’, routing = ‘‘expuh’’, e = 1, d = 200, f = range(0.5, 1.3), tau_s = range(10, 500), tau_q = range(0, 10), v_s = range(0, 1) )

    

The observed data in ‘zoo’ format IHACRES CMD version model for soil moisture accounting Exponential unit hydrograph for linear routing Fixing e and d parameter as 1 and 200 respectively Defining ranges of parameters

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

151

The Morris algorithm used to evaluate the sensitivity of NSE to parameter variation with: mm <- morris(model=evalPars, factors= names(getFreeParsRanges(model.str)), r=1000, design=list(type=’’oat’’,levels=10,grid.jump=2), binf=sapply(getFreeParsRanges(model.str),min), bsup=sapply(getFreeParsRanges(model.str),max), object=model.str, objective= hmadstat(‘‘r.squared’’)(Q, X)/(2hmadstat(‘‘r.squared’’)(Q, X))

 Utility function to evaluate parameters on a model  Names of factors/parameters  Number of repetitions  List specifying design type and its parameters  Minimum value of each non-fixed parameter  Maximum value of each non-fixed parameter  Hydromad model object to use to evaluate parameters  NSE objective function, calculated using observed (Q) and predicted (X) flows and the built-in function NSE. The tilde () indicates this is an R ‘formula’ where the values for Q and X will be substituted at the time of evaluation

)

Sensitivity of NSElog and the hybrid function (NSE + NSElog)/2 are calculated by respectively using: objective= hmadstat(‘‘r.sq.log’’)(Q, X)/(2-hmadstat(‘‘r.sq.log’’)(Q, X)) objective= 0.5⁄(hmadstat(‘‘r.squared’’)(Q, X)/(2-hmadstat(‘‘r.squared’’)(Q, X)))+0.5⁄(hmadstat(‘‘r.sq.log’’)(Q, X)/(2hmadstat(‘‘r.sq.log’’)(Q, X)))

Scalar-valued target functions that are not built in can be written as functions that refer to X (modeled flows) and Q (observed flows), e.g. the frequency of low flows (F20), thres.Q20<- as.numeric(quantile(obs$Q, probs=c(0.2), na.rm=TRUE)) F20<- function(X) length(which(X
The Sobol algorithm is used to evaluate the sensitivity of F20 to parameter variation with: n <- 10,000 X1 <- parameterSets(getFreeParsRanges(model.str),n) X2 <- parameterSets(getFreeParsRanges(model.str),n) x <- sobol2002(model = evalPars, X1 = X1, X2 = X2, nboot = 100, object=model.str, objective= F20(X) )

       

Number of random samples for parameters First random sample Second random sample Model to be analyzed Selecting parameters to use Number of bootstrap replicates Hydromad model object to use to evaluate parameters F20 prediction function

A working demonstration is also available after loading the hydromad package using the command: demo(sensitivity). References Alam, F.M., McNaught, K.R., Ringrose, T.J., 2004. Using Morris’ randomized OAT design as a factor screening method for developing simulation metamodels. In: Proceedings of the 36th Conference on Winter Simulation, Washington, DC, 2004, pp. 949–957. Alısß, Ö.F., Rabitz, H., 2001. Efficient implementation of high dimensional model representations. J. Math. Chem. 29 (2), 127–142. http://dx.doi.org/10.1023/ A:1010979129659. Anctil, F., Perrin, C., Andréassian, V., 2004. Impact of the length of observed records on the performance of ANN and of conceptual parsimonious rainfall–runoff forecasting models. Environ. Modell. Softw. 19 (4), 357–368. http://dx.doi.org/ 10.1016/S1364-8152(03)00135-X. Andréassian, V., Lerat, J., Le Moine, N., Perrin, C., 2012. Neighbors: nature’s own hydrological models. J. Hydrol. 414–415 (11), 49–58. http://dx.doi.org/10.1016/ j.jhydrol.2011.10.007. Andrews, F.T., Croke, B.F.W., Jakeman, A.J., 2011. An open software environment for hydrological model assessment and development. Environ. Modell. Softw. 26 (10), 1171–1185. http://dx.doi.org/10.1016/j.envsoft.2011.04.006. Beck, M.B., 1987. Water quality modeling: a review of the analysis of uncertainty. Water Resour. Res. 23 (8), 1393–1442. http://dx.doi.org/10.1029/ WR023i008p01393. Blasone, R.-S., Vrugt, J.A., Madsen, H., Rosbjerg, D., Robinson, B.A., Zyvoloski, G.A., 2008. Generalized likelihood uncertainty estimation (GLUE) using adaptive Markov Chain Monte Carlo sampling. Adv. Water Resour. 31 (4), 630–648. http://dx.doi.org/10.1016/j.advwatres.2007.12.003. Borgonovo, E., 2006. Measuring uncertainty importance: investigation and comparison of alternative approaches. Risk Anal. 26 (5), 1349–1361. http:// dx.doi.org/10.1111/j.1539-6924.2006.00806.x.

Burnash, R.J.C., Ferral, R.L., McGuire, R.A., 1973. A Generalized Streamflow Simulation System: Conceptual Modeling for Digital Computers. Technical Report. U.S. Natl. Weather Serv., Sacramento, Calif. Burnash, R.J.C., 1995. The NWS river forecast system: catchment modelling. In: Singh, V.P. (Ed.), Computer Models of Watershed Hydrology. Water Resour. Publ., Highlands Ranch, Colorado, pp. 311–366. Campolongo, F., Cariboni, J., Saltelli, A., 2007. An effective screening design for sensitivity analysis of large models. Environ. Modell. Softw. 22 (10), 1509–1518. http://dx.doi.org/10.1016/j.envsoft.2006.10.004. Campolongo, F., Saltelli, A., 1997. Sensitivity analysis of an environmental model: an application of different analysis methods. Reliab. Eng. Syst. Safe 57 (1), 49–69. http://dx.doi.org/10.1016/S0951-8320(97)00021-5. Cariboni, J., Gatelli, D., Liska, R., Saltelli, A., 2007. The role of sensitivity analysis in ecological modelling. Ecol. Model. 203 (1–2), 167–182. http://dx.doi.org/ 10.1016/j.ecolmodel.2005.10.045. Chiew, F.H.S., Kirono, D.G.C., Kent, D.M., Frost, A.J., Charles, S.P., Timbal, B., Nguyen, K.C., Fu, G., 2010. Comparison of runoff modelled using rainfall from different downscaling methods for historical and future climates. J. Hydrol. 387 (1–2), 10–23. http://dx.doi.org/10.1016/j.jhydrol.2010.03.025. Chiew, F.H.S., Peel, M.C., Western, A.W., 2002. Application and testing of the simple rainfall–runoff model SIMHYD. In: Singh, V.P., Frevert, D.K. (Eds.), Mathematical Models of Small Watershed Hydrology and Applications. Water Resour. Publ., Colorado, pp. 335–367. Croke, B.F.W., Jakeman, A.J., 2004. A catchment moisture deficit module for the IHACRES rainfall–runoff model. Environ. Modell. Softw. 19 (1), 1–5. http:// dx.doi.org/10.1016/j.envsoft.2003.09.001. Edijatno, Nascimento, N.O., Yang, X., Makhlouf, Z., Michel, C., 1999. GR3J: a daily watershed model with three free parameters. Hydrol. Sci. J. 44 (2), 263–277. http://dx.doi.org/10.1080/02626669909492221.

152

M.-J. Shin et al. / Journal of Hydrology 503 (2013) 135–152

Evans, J.P., Jakeman, A.J., 1998. Development of a simple, catchment-scale, rainfall– evapotranspiration–runoff model. Environ. Modell. Softw. 13 (3–4), 385–393. http://dx.doi.org/10.1016/S1364-8152(98)00043-7. Hansen, D.P., Jakeman, A.J., Kendall, C., Weizu, G., 1997. Identification of internal flow dynamics in two experimental catchments. Math. Comput. Simul. 43 (3–6), 367–375. http://dx.doi.org/10.1016/S0378-4754(97)00021-9. Homma, T., Saltelli, A., 1996. Importance measures in global sensitivity analysis of nonlinear models. Reliab. Eng. Syst. Safe 52 (1), 1–17. http://dx.doi.org/10.1016/ 0951-8320(96)00002-6. Jakeman, A.J., Hornberger, G.M., 1993. How much complexity is warranted in a rainfall–runoff model? Water Resour. Res. 29 (8), 2637–2649. http://dx.doi.org/ 10.1029/93WR00877. Jakeman, A.J., Littlewood, I.G., Whitehead, P.G., 1990. Computation of the instantaneous unit hydrograph and identifiable component flows with application to two small upland catchments. J. Hydrol. 117 (1–4), 275–300. http://dx.doi.org/10.1016/0022-1694(90)90097-H. Kim, H.S., Croke, B.F.W., Jakeman, A.J., Chiew, F.H.S., 2011. An assessment of modelling capacity to identify the impacts of climate variability on catchment hydrology. Math. Comput. Simul. 81 (7), 1419–1429. http://dx.doi.org/10.1016/ j.matcom.2010.05.007. Kucherenko, S., Tarantola, S., Annoni, P., 2012. Estimation of global sensitivity indices for models with dependent variables. Comput. Phys. Commun. 183 (4), 937–946. http://dx.doi.org/10.1016/j.cpc.2011.12.020. Le Moine, N., Andréassian, V., Mathevet, T., 2008. Confronting surface- and groundwater balances on the La Rochefoucauld-Touvre karstic system (Charente, France). Water Resour. Res. 44 (3), W03403. http://dx.doi.org/ 10.1029/2007WR005984. Littlewood, I.G., Marsh, T.J., 1996. Re-assessment of the monthly naturalized flow record for the River Thames at Kingston since 1883, and the implications for the relative severity of historical droughts. Regul. River 12 (1), 13–26. http:// dx.doi.org/10.1002/(SICI)1099-1646(199601)12:1<13::AID-RRR364>3.0.CO;2N. Mathevet, T., Michel, C., Andréassian, V., Perrin, C., 2006. A bounded version of the Nash–Sutcliffe criterion for better model assessment on large sets of basins. Large Sample Basin Experiments for Hydrological Model Parameterization: Results of the Model Parameter Experiment – MOPEX, vol. 307. IAHS Publ., pp. 211–219. Morris, M.D., 1991. Factorial sampling plans for preliminary computational experiments. Technometrics 33 (2), 161–174. http://dx.doi.org/10.1080/ 00401706.1991.10484804. Moussu, F., Oudin, L., Plagnes, V., Mangin, A., Bendjoudi, H., 2011. A multi-objective calibration framework for rainfall–discharge models applied to karst systems. J. Hydrol. 400 (3–4), 364–376. http://dx.doi.org/10.1016/j.jhydrol.2011.01.047. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I – a discussion of principles. J. Hydrol. 10 (3), 282–290. http://dx.doi.org/ 10.1016/0022-1694(70)90255-6. Norton, J.P., 2009. Selection of Morris trajectories for initial sensitivity analysis. In: 15th IFAC Symposium on System Identification (SYSID 2009). St. Malo, France. http://dx.doi.org/0.3182/20090706-3-FR-2004.00111. Nossent, J., Elsen, P., Bauwens, W., 2011. Sobol’ sensitivity analysis of a complex environmental model. Environ. Modell. Softw. 26 (12), 1515–1525. http:// dx.doi.org/10.1016/j.envsoft.2011.08.010. Peck, E.L., 1976. Catchment Modeling and Initial Parameter Estimation for the National Weather Service River Forecast System. In: NOAA Tech. Memo. NWS HYDRO-31. Hydrol. Res. Lab., Silver Spring, Md. Perrin, C., 2000. Towards an Improvement of a Lumped Rainfall–Runoff Model through a Comparative Approach. Ph.D. Thesis, Univ. Joseph Fourier, Grenoble, France (in French). Perrin, C., Michel, C., Andréassian, V., 2001. Does a large number of parameters enhance model performance? Comparative assessment of common catchment model structures on 429 catchments. J. Hydrol. 242 (3–4), 275–301. http:// dx.doi.org/10.1016/S0022-1694(00)00393-0. Perrin, C., Michel, C., Andréassian, V., 2003. Improvement of a parsimonious model for streamflow simulation. J. Hydrol. 279 (1–4), 275–289. http://dx.doi.org/ 10.1016/S0022-1694(03)00225-7. Petheram, C., Rustomji, P., Chiew, F.H.S., Vleeshouwer, J., 2012. Rainfall–runoff modelling in northern Australia: a guide to modelling strategies in the tropics. J. Hydrol. 462–463, 28–41. http://dx.doi.org/10.1016/j.jhydrol.2011.12.046. Podger, G., 2004. RRL Rainfall Runoff Library User Guide. Cooperative Research Centre for Catchment Hydrology. . Porter, J.W., 1972. The Synthesis of Continuous Streamflow from Climatic Data by Modelling with a Digital Computer. Ph.D. Thesis, Dept. of Civil Eng., Monash Univ., Clayton, Vic., Australia. Porter, J.W., McMahon, T.A., 1975. Application of a catchment model in southeastern Australia. J. Hydrol. 24 (1–2), 121–134. http://dx.doi.org/ 10.1016/0022-1694(75)90146-8. Post, D.A., Chiew, F.H.S., Teng, J., Viney, N.R., Ling, F.L.N., Harrington, G., Crosbie, R.S., Graham, B., Marvanek, S., McLoughlin, R., 2012. A robust methodology for

conducting large-scale assessments of current and future water availability and use: a case study in Tasmania, Australia. J. Hydrol. 412–413, 233–245. http:// dx.doi.org/10.1016/j.jhydrol.2011.02.011. Pujol, G., Iooss, B., Janon, A., 2012. Package ‘Sensitivity’: Sensitivity Analysis. R Package Version 1.5. . Pushpalatha, R., Perrin, C., Le Moine, N., Andréassian, V., 2012. A review of efficiency criteria suitable for evaluating low-flow simulations. J. Hydrol. 420–421, 171– 182. http://dx.doi.org/10.1016/j.jhydrol.2011.11.055. Saltelli, A., 2002. Making best use of model evaluations to compute sensitivity indices. Comput. Phys. Commun. 145 (2), 280–297. http://dx.doi.org/10.1016/ S0010-4655(02)00280-1. Saltelli, A., Annoni, P., 2010. How to avoid a perfunctory sensitivity analysis. Environ. Modell. Softw. 25 (12), 1508–1517. http://dx.doi.org/10.1016/ j.envsoft.2010.04.012. Saltelli, A., Chan, K., Scott, E.M., 2000. Sensitivity Analysis. Wiley, New York. Silberstein, R.P., Aryal, S.K., Durrant, J., Pearcey, M., Braccia, M., Charles, S.P., Boniecka, L., Hodgson, G.A., Bari, M.A., Viney, N.R., McFarlane, D.J., 2012. Climate change and runoff in south-western Australia. J. Hydrol. 475, 441–455. http:// dx.doi.org/10.1016/j.jhydrol.2012.02.009. Sobol, I.M., 1990. Sensitivity estimates for nonlinear mathematical models. Matem. Modelirovanie 2 (1), 112–118 (in Russian). Sobol, I.M., 1993. Sensitivity analysis for nonlinear mathematical models. Math. Model. Comput. Exp. 1 (4), 407–414. Sobol, I.M., 2001. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math. Comput. Simulat. 55 (1–3), 271–280. http:// dx.doi.org/10.1016/S0378-4754(00), 00270-6. Sun, X.Y., Newham, L.T.H., Croke, B.F.W., Norton, J.P., 2012. Three complementary methods for sensitivity analysis of a water quality model. Environ. Modell. Softw. 37, 19–29. http://dx.doi.org/10.1016/j.envsoft.2012.04.010. Thyer, M., Renard, B., Kavetski, D., Kuczera, G., Franks, S.W., Srikanthan, S., 2009. Critical evaluation of parameter consistency and predictive uncertainty in hydrological modeling: a case study using Bayesian total error analysis. Water Resour. Res. 45 (12), W00B14. http://dx.doi.org/10.1029/2008WR006825. van Griensven, A., Meixner, T., Grunwald, S., Bishop, T., Diluzio, M., Srinivasan, R., 2006. A global sensitivity analysis tool for the parameters of multi-variable catchment models. J. Hydrol. 324 (1–4), 10–23. http://dx.doi.org/10.1016/ j.jhydrol.2005.09.008. van Werkhoven, K., Wagener, T., Reed, P., Tang, Y., 2008. Characterization of watershed model behavior across a hydroclimatic gradient. Water Resour. Res. 44 (1), W01429. http://dx.doi.org/10.1029/2007WR006271. van Werkhoven, K., Wagener, T., Reed, P., Tang, Y., 2009. Sensitivity-guided reduction of parametric dimensionality for multi-objective calibration of watershed models. Adv. Water Resour. 32 (6), 1154–1169. http://dx.doi.org/ 10.1016/j.advwatres.2009.03.002. Vaze, J., Post, D.A., Chiew, F.H.S., Perraud, J.M., Viney, N.R., Teng, J., 2010. Climate non-stationarity – validity of calibrated rainfall–runoff models for use in climate change studies. J. Hydrol. 394 (3–4), 447–457. http://dx.doi.org/ 10.1016/j.jhydrol.2010.09.018. Viney, N.R., Bormann, H., Breuer, L., Bronstert, A., Croke, B.F.W., Frede, H., Gräff, T., Hubrechts, L., Huisman, J.A., Jakeman, A.J., Kite, G.W., Lanini, J., Leavesley, G., Lettenmaier, D.P., Lindström, G., Seibert, J., Sivapalan, M., Willems, P., 2009. Assessing the impact of land use change on hydrology by ensemble modelling (LUCHEM) II: ensemble combinations and predictions. Adv. Water Resour. 32 (2), 147–158. http://dx.doi.org/10.1016/j.advwatres.2008.05.006. Wang, J., Li, X., Lu, L., Fang, F., 2013. Parameter sensitivity analysis of crop growth models based on the extended Fourier Amplitude Sensitivity Test method. Environ. Modell. Softw. 48, 171–182. http://dx.doi.org/10.1016/ j.envsoft.2013.06.007. Wagener, T., McIntyre, N., Lees, M.J., Wheater, H.S., Gupta, H.V., 2003. Towards reduced uncertainty in conceptual rainfall runoff modelling: dynamic identifiability analysis. Hydrol. Process. 17 (2), 455–476. http://dx.doi.org/ 10.1002/hyp.1135. Yang, J., 2011. Convergence and uncertainty analyses in Monte-Carlo based sensitivity analysis. Environ. Modell. Softw. 26 (4), 444–457. http://dx.doi.org/ 10.1016/j.envsoft.2010.10.007. Yapo, P.O., Gupta, H.V., Sorooshian, S., 1996. Automatic calibration of conceptual rainfall–runoff models: sensitivity to calibration data. J. Hydrol. 181 (1–4), 23– 48. http://dx.doi.org/10.1016/0022-1694(95)02918-4. Yatheendradas, S., Wagener, T., Gupta, H., Unkrich, C., Goodrich, D., Schaffner, M., Stewart, A., 2008. Understanding uncertainty in distributed flash flood forecasting for semiarid regions. Water Resour. Res. 44, W05S19, http:// dx.doi.org/10.1029/2007WR005940. Ye, W., Bates, B.C., Viney, N.R., Sivapalan, M., Jakeman, A.J., 1997. Performance of conceptual rainfall–runoff models in low-yielding ephemeral catchments. Water Resour. Res. 33 (1), 153–166. http://dx.doi.org/10.1029/96WR02840. Zhang, Y., Chiew, F.H.S., 2009. Relative merits of different methods for runoff predictions in ungauged catchments. Water Resour. Res. 45 (7), W07412. http:// dx.doi.org/10.1029/2008WR007504.