Comparing methods for estimating flow duration curves at ungauged sites

Journal of Hydrology 434–435 (2012) 78–94

Contents lists available at SciVerse ScienceDirect

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

Comparing methods for estimating flow duration curves at ungauged sites D.J. Booker ⇑, T.H. Snelder National Institute of Water and Atmospheric Research, PO Box 8602, Riccarton, Christchurch, New Zealand

a r t i c l e

i n f o

Article history: Received 26 July 2011 Received in revised form 8 January 2012 Accepted 14 February 2012 Available online 22 February 2012 This manuscript was handled by Andras Bardossy, Editor-in-Chief, with the assistance of Erwin Zehe, Associate Editor Keywords: Flow duration curves Ungauged sites Probability distribution functions Random forests Mixed-effects Dissimilarity modelling

s u m m a r y Flow duration curves (FDCs) are a useful tool for characterising hydrological regimes and flow variability. FDCs observed at 379 gauging stations located across New Zealand were analysed with the aim of investigate how parameterisation and generalisation combine to influence the accuracy of empirically predicted FDCs at ungauged sites. The appropriateness of four strategies for estimating FDCs was compared: (a) parameterise then generalise; (b) parameterise then regionalise then generalise; (c) parameterise and generalise together; and (d) FDC substitution. These strategies were deployed using various combinations of methods for calculating parameters that describe the shape of FDCs (polynomial expressions and probability distribution functions) and then methods for estimating these parameters at ungauged sites using available catchment characteristics (stepwise linear regression and random forests). A parameterise and generalise together strategy was devised by applying a mixed-effects approach. A jack-knife cross-validation procedure was used to provide an independent test of each method for estimating the FDC at ungauged sites. For parameterise then regionalise strategies, it was found that the combination of parameterisation method and generalisation method together, rather than either in isolation, was important in determining overall performance. Results indicated that predictive capability varied between methods and across exceedence percentiles. The mixed-effects approach provided the most parsimonious method for estimating FDC at ungauged sites. A method using the generalised extreme value probability distribution that was generalised using random forests was the most accurate method of estimating flow duration curves at ungauged sites across New Zealand. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction The flow duration curve (FDC) is a tool used to describe hydrological regimes. The FDC represents the relationship between magnitude and frequency of flow by defining the proportion of time for which any discharge is equalled or exceeded (Vogel and Fennessey, 1994). FDCs represent useful graphical and analytical tools for evaluating flow variability at a particular site. Information contained within the FDC can be used for water resource assessments including hydropower design schemes (Warnick, 1984), reliability of water supply (McMahon, 1993), water quality assessments (Vogel and Fennessey, 1995) and the evaluation of river habitats (Booker and Dunbar, 2004). FDCs can be calculated from historical discharge time-series’ where records of sufficient length are available. However, these records are only available where the required data have been observed at gauging stations, and water resource assessments are often required for ungauged locations (Nathan and McMahon, 1992). This creates a need for estimation of hydrological statistics including

⇑ Corresponding author. Tel.: +64 (0)3 348 8987; fax: +64 (0)3 348 7891. E-mail addresses: [email protected] (D.J. Booker), [email protected] (T.H. Snelder). 0022-1694/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2012.02.031

the FDC at ungauged sites. Examples of attempts to estimate FDCs at ungauged sites include those applied in Canada (LeBoutillier and Waylen, 1993), France (Sauquet and Catalogne, 2011), Greece (Mimikou and Kaemaki, 1985), India (Singh et al., 2001), Italy (Franchini and Suppo, 1996; Castellarin et al., 2004), Turkey (Cigizoglu and Bayazit, 2000), Philippines (Quimpo et al., 1983), Portugal (Croker et al., 2003), South Africa (Smakhtin et al., 1997), Switzerland (Ganora et al., 2009), Taiwan (Yu et al., 2002), United Kingdom (Holmes et al., 2002) and United States (Fennessey and Vogel, 1990). Many of these empirical approaches comprise two steps. First, the shapes of many FDCs are described by calculating values of statistical parameters, a process referred to herein as parameterisation. Second, the shape parameters are related to catchment characteristics. This allows estimation of FDCs at ungauged sites. In this paper we refer to this process as generalisation. Therefore, empirical methods for estimation of FDCs at ungauged sites typically comprise parameterisation and generalisation. Various methods have been used to estimate FDCs for ungauged sites, including: regional regression approaches (e.g., Nathan and McMahon, 1992); defining regional prediction curves (e.g., Smakhtin et al., 1997); mapping and interpolation of dimensionless flow indices (e.g., Vandewiele and Elias, 1995; Arnell, 1995); deriving properties of synthetic time-series (e.g., Clausen et al., 1994); and

79

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

empirical orthogonal functions in combination with a hydrological regionalisation (Sauquet and Catalogne, 2011). Smakhtin (2001) provides a detailed review of these procedures. FDCs for ungauged sites can also be estimated by substituting the unknown FDC with an available observed FDC. For example, Ganora et al. (2009) used a dissimilarity method whereby a matrix describing dissimilarities between all pairs of FDCs was related to a matrix describing dissimilarities between the catchment characteristics of all pairs of gauging stations. This study focused on empirical methods for estimating FDCs at ungauged sites, rather than process-based methods that seek to transfer parameters of physically-based hydrological models (e.g., Wagener and Wheater, 2006; Bárdossy, 2007; Buytaert and Bevan, 2009). The overall aim was to investigate how parameterisation and generalisation combine to influence the accuracy of empirically predicted FDCs at ungauged sites. In order to fulfil this aim the appropriateness of several empirical methods for estimating FDCs at ungauged sites were investigated. The specific objectives of the study were to: (a) compare the accuracy and parsimony of several methods of statistical description of FDCs observed at many gauging stations across New Zealand; (b) compare methods for parameterising these statistical descriptions from catchment characteristics for application at ungauged sites; (c) identify the most accurate method of FDC-substitution for application at ungauged sites; and (d) quantify the uncertainty associated with estimated FDCs for ungauged sites.

100 km N

O

40 S

Cool Dry Cool Wet Cool Extremely Wet Warm Dry Warm Wet Warm Extremely Wet

2. Data description 2.1. Hydrological data and flow duration curves FDCs were calculated using mean daily flows observed at 379 gauging stations with publicly available records of five full years or longer. Only gauging stations whose catchments were natural or had only minimal abstraction and impoundment were included in the analysis. See Snelder et al. (2005a) for further details on gauging station selection. These gauging stations were located throughout New Zealand (Fig. 1), and represented a wide range of hydrological conditions (Table 1, Table 2). The observed timeseries’ did not all cover the same time periods (Fig. 2). To allow comparisons between the shapes of FDCs, each was standardised by dividing by the mean flow for that gauging station. 2.2. Catchment characteristics A GIS representation of the New Zealand river network comprising 550,000 segments, their unique upstream catchments and an associated database of catchment characteristics were used to provide information for each gauging station. The catchment characteristics include a range of categorical and continuous variables (e.g., Table 1) (Snelder and Biggs, 2002; Snelder et al., 2004; Leathwick et al., 2011). The GIS river network and associated databases have previously been used to define a hierarchical classification of New Zealand’s rivers called the River Environment Classification (REC; Snelder and Biggs, 2002). These databases provide inventories for river resource analysis and management purposes (Snelder and Hughey, 2005; Leathwick et al., 2011; Clapcott et al., 2010, 2011). They have also been used to create nationwide models for estimating flow statistics such as flood flows (Pearson and McKerchar, 1989), low flows (Pearson, 1995) and mean flow (Woods et al., 2006) at ungauged sites using relationships between these hydrological metrics and catchment characteristics. Snelder et al. (2005b) showed that grouping river segments by nested categorical subdivisions of climate and topography, known as the Sourceof-Flow grouping factor (Table 1), provided an a priori hydrological regionalisation.

O

175 E

Fig. 1. Map showing the locations and climate category of the gauging stations used in this study.

3. Strategies Four strategies were employed for estimating flow duration curves at ungauged sites: (a) parameterise then generalise (PG); (b) parameterise then regionalise then generalise (PRG); (c) parameterise and generalise together (PGT); and (d) FDC substitution. Variations in the statistical approaches within these strategies lead to the definition of 19 methods which are fully described below and summarised in Table 3. 3.1. Parameterise then generalise PG strategies comprise two separate steps. In the first step predefined characteristics of each observed FDC are quantified by calculating parameter values using parametric methods. This step can be accomplished using regression methods (Fennessey and Vogel, 1990) or by calculating linear moments (Hosking, 1990). In the second step the values calculated for each parameter for all observed FDCs are related to catchment characteristics. This provides a method for estimating FDCs at ungauged sites from known catchment characteristics. 3.2. Parameterise then regionalise then generalise Hydrological regionalisations attempt to classify river network locations into groups with similar hydrological characteristics (Burn, 1997). Regionalisations can be created either by defining regions a priori (e.g., Snelder and Biggs, 2002) or by classifying

80

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Table 1 Summary of the defining characteristics, categories and category membership criteria that combine to define Source-of-Flow groupings within the REC.

a b

Defining characteristic

Categories

Notation

Category membership criteria

Climate

Warm-extremely-wet Warm-wet Warm-dry Cool-extremely-wet Cool-wet Cool-dry

WX WW WD CX CW CD

Warm: mean annual temperature P 12 °C Cool: mean annual temperature < 12°C Extremely Wet: mean annual effective precipitationa P 1500 mm Wet: mean annual effective precipitation > 500 and < 1500 mm Dry: mean annual effective precipitation 6 500 mm

Topography

Glacial-mountain Mountain Hill Low-elevation Lake

GM M H L Lk

GM: M and % permanent ice > 1.5% M: > 50% annual rainfall volume above 1000 m H: 50% rainfall volume between 400 and 1000 m L: 50% rainfall below 400 m Lk: Lake influence indexb>0.033

Effective precipitation = annual rainfall  annual potential evapotranspiration. See Snelder and Biggs (2002) for a description.

Table 2 Source-of-Flow groupings within the REC for gauging stations used in this study. Climate

Topography Glacialmountain

Warm Extremely Wet Warm Wet Warm Dry Cool Extremely Wet Cool Wet Cool Dry

Mountain

Hill

Lowland

catchment characteristics and then finding the nearest gauging station to the ungauged site; and dissimilarity modelling (e.g., Ganora et al., 2009).

Lake

0

0

1

2

0

0 0 8 2 0

0 0 14 28 1

0 0 39 73 40

68 12 6 31 36

0 0 13 4 1

natural flow regimes according to their similarity with respect to hydrological indices (e.g., Isik and Singh, 2008; Snelder et al., 2009; Sadri and Burn, 2011). PRG strategies are those that apply separate parameterise then generalise procedures within groups of a pre-defined hydrological regionalisation. 3.3. Parameterise and generalise together One difficulty with PG strategies, regardless of regionalising, arises because the parameters used to describe the FDCs are not necessarily independent, but they may be generalised independently. This can lead to combinations of parameter values being calculated for new sites that, when combined, produce spurious FDC shapes. The same phenomena may occur in physically-based hydrological modelling where interpolation of individual parameter values can lead to unreasonable model parameters and results (Bárdossy, 2007). One method for dealing with uncertainties in the parameter transformation process is to relate parameter behaviour to catchment properties and intercatchment similarities by creating model ensembles (Buytaert and Bevan, 2009). In the case of empirical estimation of FDCs, independent generalisation of parameter values could produce FDCs that are not monotonically increasing. One method of overcoming this difficulty is to employ a PGT strategy, such as a mixed-effects approach (Pinheiro and Bates, 2000) in which all FDCs are parameterised and generalised together within the same procedure though specification of both the fixed and random effects. 3.4. FDC substitution An alternative approach to estimating FDCs for ungauged sites is to substitute with an observed FDC. Various methods can be used to identify which observed FDC should be used as the substitute. These methods include: finding the nearest gauging station to the ungauged site; isolating gauges that might have similar

4. Methods 4.1. Parameterise then generalise Two different statistical approaches to parameterising an FDC were compared for the first step in PG strategies. These approaches were: (a) fitting of various linear regression models; and (b) calculating linear moments and then generating an FDC using various probability distribution functions. For standardised flows the first linear moment ,l1, represents the mean and is always equal to 1, the second moment , l2, represents the variance, the third moment, lca, represents skewness and the fourth moment, lkur, represents kurtosis (Hosking, 1990; Vogel and Fennessey, 1993). For the second step in PG strategies, two statistical approaches for generalisation were compared. These approaches were: (a) stepwise linear regression; (b) random forests.

4.1.1. Linear regression models Several polynomial equations were fitted to the observed FDCs. For each of the 379 observed FDCs, j, discharge, Q, was standardised by dividing by the mean and then logged to the base 10. Log standardised Q was modelled as a function of exceedance percentile by applying five linear regression models each with a different formulation (Table 4). For these linear regressions percent of time flow was not exceeded was transformed to be normal-reduced, Zi (Fennessey and Vogel, 1990). Use of Zi  Zmin (where Zmin was a constant of 4.66) allowed a base 10 log transformation of the already normal-reduced exceedance percentiles. Use of Zi  1 allowed the FDC to be centred within the range of percentiles, and therefore the possibility for odd polynomial terms to describe asymptotes at both high and low flows. A Zi  1 transformation centres the FDC around the 16th flow percentile. In further discussion Eqs. (1)–(5) (Table 4) are referred to as the Linear, Log, Cubic, Squared and Comb3 models respectively. Each linear model was fit using ordinary least squares linear regression with the assumptions of normal errors and constant variance (Chambers, 1992). The different linear models were compared to determine which of Eqs. (1)–(5) provided the most accurate, yet parsimonious model of the observed FDCs. For each site and each model the Akaike information criterion (AIC), residual standard error and adjusted r-squared were calculated and compared. AIC is a measure of the trade-off between degrees of freedom and fit of the model (Akaike, 1973).

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

350

terms on the basis of AIC alone has been shown to be somewhat liberal in its choice of terms, a value of k = 4 was used for the multiple of the number of degrees of freedom used for the penalty in the stepwise procedure (Venables and Ripley, 2002). Minimal adequate models obtained from stepwise reduction of Eq. (6) for each observed FDC are referred to as Step models in further discussion.

200

250

300

4.1.2. Probability distribution functions Seven probability distribution functions were identified as possible methods for generating FDCs as defined by the distribution of standardised discharge. These probability distribution functions were based on various pre-defined probability distribution functions and their transformations (Table 5), which are themselves calculated from linear moments (Laio et al., 2009). For each of the 379 observed FDCs, the AIC and the Anderson Darling Criteria (ADC) were calculated for each probability distribution following the method of Laio et al. (2009). These two different model selection criteria were both designed to assess the appropriateness of a distribution to represent hydrological frequency data. The distribution of both AIC and ADC calculated for each probability distribution and for each observed FDC were compared. This comparison was made for all FDCs together and after having grouped FDCs by Source-of-Flow in order to assess the appropriateness of each probability distribution function to summarise FDC characteristics across a range of hydrological settings. The majority of the 379 observed FDCs exhibited permanent flow, but 33 contained zero flows. The presence of zero flows can be problematic for probability distribution functions requiring log transformations. In these cases zero flows were set to be one tenth of lowest observed non-zero flow for each FDC.

50

100

150

Sites

81

1950

1960

1970

1980

1990

2000

2010

Year Fig. 2. Time-periods for which data were available from each gauging station. Note seven records started before 1950.

As an additional aid to finding the most appropriate nationalscale equation from the candidate equations, a saturated model was also defined that included many terms that could possibly combine to best model the FDC:

logðQ i =Q Þ ¼ a0 þ a1 ðZ i  1Þ þ a2 logðZ i  Z min Þ þ a3 ðZ i  Z min Þ2 þ a4 ðZ i  1Þ3

ð6Þ

Standard forwards and backwards stepwise linear regression was applied to Eq. (6) to identify the minimal adequate model from the terms included in this saturated model for each FDC separately. The Akaike information criterion (AIC; Akaike, 1973) was used to apply a penalised log likelihood method to evaluate the trade-off between degrees of freedom and fit of the model as explanatory parameters are added or removed (Crawley, 2002). As selecting

4.1.3. Generalisation Each of the parameters from each of the various linear regression models (Eqs. (1)–(5); Table 4), along with the linear moments of raw flow data (required for Eqs. (7)–(10); Table 5) and linear moments of the log transformed flow data (required for Eqs. (11)–(13); Table 5) was modelled independently as a function of available catchment characteristics. This provided several methods for estimating FDCs at ungauged sites. For each parameter the distribution of values fitted to the observed FDCs was modelled as a function of a suite of available independent continuous variables describing the physical characteristics of the upstream catchment (Table 6). These variables were chosen to include characteristics that were likely to influence hydrological processes, but no attempt was made to explicitly model any hydrological processes. Two different statistical methods were used to model each fitted parameter separately. To avoid any subjectivity in the fitting process we choose to use only fully automated statistical techniques. The first generalisation method used stepwise multiple-linear regression. This method was used to identify the minimal adequate linear model of each parameter as a function of the candidate explanatory variables including all two-way interactions (Table 6). As above, AIC with a value of k = 4 was used to apply a penalised log likelihood method to evaluate the trade-off between degrees of freedom and fit of the model as more explanatory parameters are added into it (Crawley, 2002; Venables and Ripley, 2002). Completion of this process created a statistical model that could be used for estimating each parameter at an ungauged location. The second generalisation method used random forests (Breiman, 2001). This method uses machine-learning to combine many regression trees to produce more accurate regressions (Cutler et al., 2007). Random forests were used to model each parameter as a function of the explanatory variables (Table 6). A Random Forest model comprises an ensemble of regression trees (a forest) from which a final prediction is based on the predictions averaged over

82

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Table 3 Descriptions of methods used to estimate FDCs. Strategy

Method name

Description

Parameterise then generalise

Linear random forest

Parameterise using a straight line regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded. Generalisation using random forests Parameterise using a straight line regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded. Generalisation using stepwise reduction of linear models with two-way interactions Parameterise using a polynomial regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded, which includes a quadratic and log term. Generalisation using random forests Parameterise using a polynomial regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded, which includes a quadratic and log term. Generalisation using stepwise reduction of linear models with two-way interactions Parameterise by calculating linear moments of standardised flow. Generalise using random forests. Generate FDC using the GEV probability distribution function Parameterise by calculating linear moments of standardised flow. Generalisation using stepwise reduction of linear models with two-way interactions. Generate FDC using the GEV probability distribution function Parameterise by calculating linear moments of log standardised flow. Generalise using random forests. Generate FDC using the LP3 probability distribution function Parameterise by calculating linear moments of log standardised flow. Generalisation using stepwise reduction of linear models with two-way interactions. Generate FDC using the LP3 probability distribution function

Linear step

Comb3 RandomForest

Comb3 Step

GEV RandomForest GEV Step

LP3 RandomForest LP3 Step

Parameterise then regionalise then generalise Linear Reg Then Gen

Comb3 Reg Then Gen

GEV Reg Then Gen

LP3 Reg Then Gen

Parameterise using a straight line regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded. Regionalise using Source-of-Flow classification. Generalisation by calculating the average in each class Parameterise using a polynomial regression of log standardised flow against normal-reduced percentage of time that flow is not exceeded, which includes a quadratic and log term. Regionalise using Source-of-Flow classification. Generalisation by calculating the average in each class Parameterise by calculating linear moments of standardised flow. Regionalise using Source-ofFlow classification of the REC. Generalisation by calculating the average in each class. Generate FDC using the GEV probability distribution function. Parameterise by calculating linear moments of log standardised flow. Regionalise using Sourceof-Flow classification. Generalisation by calculating the average in each class. Generate FDC using the LP3 probability distribution function

Parameterise and generalise together

Mixed

Apply a mixed-effects model. Fixed effects are a polynomial of log standardised flow against normal-reduced percentage of time that flow is not exceeded, which includes a quadratic and log term. Random effects are stream order within Source-of-Flow category

FDC substitution

Nearest REC Distance matrix Dissimilarity

Find the nearest site (Euclidean distances) within the same Source-of-Flow class Substitute using the FDC with the least dissimilar FDC shape parameters. Substitute using the FDC that a dissimilarity model predicts to be the least dissimilar, where dissimilarities between in FDC shape parameters are estimated from dissimilarities between catchment characteristics Substitute using a randomly selected FDC

Random gauge

Table 4 Definitions of equations used to model flow duration curves. Where Z is normalreduced percent of time flow is exceeded, Q is flow and c is 1. Acronym

Linear function

Equation nos.

Linear

logðQ i =Q Þ ¼ a0 þ a1 ðZ i  cÞ

(1)

Log

logðQ i =Q Þ ¼ a0 þ a1 logðZ i  Z min Þ

(2)

Cubic

logðQ i =Q Þ ¼ a0 þ a1 ðZ i  cÞ þ a2 ðZ i  cÞ3

(3)

Squared

logðQ i =Q Þ ¼ a0 þ a1 ðZ i  cÞ þ a2 ðZ i  Z min Þ2

(4)

Comb3

logðQ i =Q Þ ¼ a0 þ a1 ðZ i  cÞ

(5)

þa2 ðZ i  Z min Þ2 þ a3 ðZ i  cÞ3

all trees (Breiman, 2001; Cutler et al., 2007). A random forest model is created by drawing several bootstrap samples from the original training data and fitting a single classification tree to each sample. Independent predictions (i.e. independent of the model fitting procedure) are made for each tree from the observations that were excluded from the bootstrap sample (the OOB samples). These predictions are aggregated over all trees (the OOB predictions) and provide an estimate of the predictive performance of the model for new cases (Breiman, 2001). By-products of the random forest calculations include measures of variable importance, which are evaluated by randomly permuting each predictor variable in turn and predicting the response for the OOB observations. The decrease in prediction performance is the measure of

Table 5 Definitions of probability distribution functions used to model flow duration curves (after Laio et al. (2009)). Distribution description

Acronym

Cumulative distribution function (G) or probability distribution function (g)

Equation nos.

Gumbel or extreme value type I Normal or Gaussian

GUMBEL NORM

G(x, #) = exp[-exp(-(x - #1)/#2)] h i pffiffiffiffiffiffiffi gðx; #Þ ¼ ð1= 2p#2 Þ exp 1=2ððx  #1 Þ=#2 Þ2

(7) (8)

Generalised extreme value

GEV

Gðx; #Þ ¼ exp½ð1  ð#3 ðx  #1 ÞÞ=#2 Þ1=#3 

(9)

Gamma or Pearson type III

P3

(10)

Frechet or log transformed Gumbel Log transformed NORM Log transformed P3

EV2 LN LP3

gðx; #Þ ¼ ½1=ðj#2 jCð#3 ÞÞððx  #1 Þ=#2 Þð#3 1Þ expð½ðx  #1 Þ=#2 Þ Log transformation of Eq. (7) Log transformation of Eq. (8) Log transformation of Eq. (10)

(11) (12) (13)

83

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94 Table 6 Variables used to generalise the shape of flow duration curves for prediction at ungauged locations. Variable name

Description

LogusAnRainVar

Coefficient of variation of annual catchment rainfall (m) (logged) Mean minimum July air temperature (°C) Mean January air temperature (°C) Catchment rain days, greater than 10 mm/month (days/ year) (logged) Catchment rain days, greater than 50 mm/month (days/ year) (logged) Annual potential evapotranspiration of catchment (mm) Average elevation in the upstream catchment (m) (logged) Lake index (dimensionless) Catchment average of calcium (ordinal scale) (logged) Catchment average of hardness, induration (ordinal scale) Catchment average of particle size (ordinal scale) Catchment average of slope (m/m) Catchment average of winter solar radiation (W/m2) Catchment average of summer solar radiation (W/m2) Catchment area (m2) (logged) Stream order (Strahler stream order)

usAvTCold usAvTWarm LogusRainDays10 LogusRainDays50 usPET LogusCatElev usLake LogusCalc usHard usParticleSize usAveSlope usSolarRadWin usSolarRadSum LogCATCHAREA ORDER

questions to be answered than PG strategies, which do not consider the nested nature of the data (Gelman and Hill, 2007; Booker and Dunbar, 2008; Booker, 2010). In this context, the results take the form of fixed effects and random effects. Fixed effects model the overall response of Q in relation to Z. Random effects model variation around the overall response. Mixed-effects models were formulated using Eq. (1) to specify fixed effects, and allowed variations in each of the coefficients of Eq. (1) to vary within a hierarchical grouping structure:

Level 1ðpercentileÞ i logðQ ijkl =Q jkl Þ ¼ a0 þ a1 ðZ i  cÞ þ eikl Level 2 (Source-of-Flow) k

intercept

a0kl ¼ a0l þ U 0kl

slope

a1kl ¼ a1l þ U 1kl

intercept random component U 0kl  Nð0; s20 Þ slope random component

U1kl  Nð0; s21 Þ

Level 3 (stream order) l

intercept

a0l ¼ a0 þ U 0l

slope

a1l ¼ a1 þ U 1l

intercept random component U 0l  Nð0; s20 Þ slope random component

U 1l  Nð0; s21 Þ

Covariances and residual importance of the original variable. Importance represents the contribution to accuracy of independent predictions for each explanatory variable and is equivalent to the error resulting from dropping a term from a linear model. Each random forest was developed by growing 500 trees. As the number of trees (k) increases the generalisation error always converges, it was assumed that 500 was sufficiently high to ensure convergence. 4.2. Parameterise then regionalise then generalise The PG method could have been applied separately within hydrological regions as pre-defined by a modified version of the Source-of-Flow classes (Table 1). However, running regressions (using either linear regression or random forests) with less data points than independent variables is nonsensical. Therefore, stations were grouped into classes and the mean of each required FDC shape parameter calculated for each group. These group-mean parameters were then used to generate a representative FDC for each group. Stations were grouped into hydrological regions as defined by the Source-of-Flow classes (Table 2). For this procedure the Source-of-Flow classes were modified such that no group contained less than five stations. To do this: one Cool Dry Lake was placed with the Cool Wet Lakes; one Cool Dry Mountain was placed with the Cool Wet Mountains; two Cool Wet Glacials were placed with the Cool Wet Mountains; and one Warm Extremely Wet Hill and two Warm Extremely Wet Lowlands were placed with the Warm Wet Lowlands. 4.3. Parameterise and generalise together Linear multilevel models (Snijders and Bosker, 1999) also known as mixed-effects models (Pinheiro and Bates, 2000) were applied to parameterise a linear formulation describing observed FDCs. Mixed-effects models allow consideration of a nested hierarchy of variation in the response variable, and explanatory variables at various levels of that hierarchy. Multilevel models account for correlation of data within each level of the hierarchy. Correct degrees of freedom are retained and unbalanced data (e.g., unequal numbers of sites within classes) do not bias results. Although the results of multilevel models are initially more difficult to interpret, they provide more information and allow more

covariance covariance residual

cov ðU 0kl ; U 1kl Þ ¼ s cov ðU 0l ; U 1l Þ ¼ m jklm  nNð0; r2 Þ

ð14Þ

Several different hierarchical grouping structures could have been used to describe between-station variability within the FDC data, but we used a grouping structure which nested stream order, k, within the REC Source-of-Flow grouping factor, l (Table 1). It was hypothesised that this grouping structure was likely to represent differences in FDC shapes across the landscape, because Sourceof-Flow incorporates both climatic and topographic information and because hydrological differences are likely between catchments of different sizes (e.g., Wiltshire, 1986; Burn, 1997). This model was compared with several reduced models by removing various random effects. The formulation of the random effects was assessed by comparing AIC (Akaike, 1973) calculated for each model. P-values for Wald tests on individual (fixed effects) parameters were also examined. 4.4. FDC substitution 4.4.1. Nearest in REC class For this method, the unknown FDC at an ungauged site was substituted with an observed FDC by identifying the nearest gauging station belonging to the same Source-of-Flow grouping. The nearest gauging station was defined from Euclidean distances calculated from geographic coordinates. This method therefore used information on mapped locations of sites in combination with an a priori regionalisation of sites based on their catchment characteristics (i.e. Table 1). 4.4.2. Dissimilarity modelling Dissimilarity modelling involves the regression of a response matrix containing dissimilarities between FDCs calculated from all gauging stations, and several explanatory matrices containing the equivalent dissimilarities between gauging stations based on catchment characteristics (Legendre et al., 1994). Specifically, for n gauging stations there are n(n  1)/2 dissimilarities that correspond to the off-diagonal entries in a dissimilarity matrix. We defined dissimilarities between all pairs of Linear moments (l2, lca and lkur) calculated from observed FDCs using the Manhattan distance

84

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

metric. Although more complex regression methods can be used (e.g., Ferrier et al., 2007; Lichstein, 2007), we used multiple linear regression to relate the FDC dissimilarities to the dissimilarities calculated from the catchment characteristics (Table 6). Thus, the model expressed the FDC dissimilarity (l) as a function of the catchment characteristic dissimilarities in the form;

l ¼ b1 þ b1 g1 þ b2 g2 þ b3 g3    þ bn gn

ð15Þ

where gi are dissimilarities derived from the ith catchment characteristic. An improvement to the dissimilarity modelling method of (Ganora et al., 2009) was included to account for non-linear relationships between FDC dissimilarity and dissimilarities in catchment characteristics. The formulation shown in Eq. (15) models dissimilarities between FDCs as linear combinations of dissimilarities between catchment characteristics. However, the rate of change in dissimilarity between gauging stations is likely to vary in a non-linear manner with respect to changes in catchment characteristics. For example, dissimilarity between FDCs is likely to increase rapidly with change in catchment size for small catchments but this rate of change is likely to decrease as catchments become larger. The dissimilarity model fitting procedure tested a range of parametric transformations of the catchment characteristic dissimilarities and chose those with the strongest correlation with the response. These transformations consisted of applying a log to the base 10, applying a double log transformation and expansion of the original values either by raising them to the power of 2 or 3 (Snelder et al., 2009). Model fitting used a forward stepwise procedure to iteratively add catchment characteristics to the model. A five-fold cross validation (CV) was used to determine model complexity, i.e., to deter-

Table 7 Percentage of terms maintained in the model by the stepwise procedure by Source-ofFlow. The intercept, a0, was retained in all cases. Source-of-Flow grouping

Parameter a1(Z  1)

a2log(Z  Zmin)

a3(Z  Zmin)2

a4(Z  1)3

Cool Dry Hill Cool Dry Lake Cool Dry Lowland Cool Dry Mountain Cool Extremely Wet Glacial Cool Extremely Wet Hill Cool Extremely Wet Lake Cool Extremely Wet Lowland Cool Extremely Wet Mountain Cool Wet Glacial Cool Wet Hill Cool Wet Lake Cool Wet Lowland Cool Wet Mountain Warm Dry Lowland Warm Extremely Wet Hill Warm Extremely Wet Lowland Warm Wet Lowland

88 100 94 100 50

92 100 97 100 88

92 100 86 100 100

85 100 94 100 100

92

95

100

100

100

77

92

100

100

100

100

100

93

79

93

100

100 90 100 94 86 100 100

100 93 100 90 96 75 100

100 89 100 87 96 92 100

100 93 50 100 82 100 100

50

100

100

100

87

90

84

87

mine the justifiable number of predictors to include in the model and prevent over fitting (Snelder et al., 2009). In this CV procedure, five mutually exclusive subsets of the full dataset each containing

Stepwise

Model name

Comb3 Squared Cubic Log Linear -1000

-500

0

AIC

Stepwise

Model name

Comb3 Squared Cubic Log Linear 0.70

0.75

0.80

0.85

Adjusted r

0.90

0.95

1.00

2

Fig. 3. Distributions of AIC and adjusted r-squared for all sites (n = 379) and various model formulations.

85

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Warm Extremely Wet_Lowland

Warm Wet_Lowland

PDF

Order

Generalized extreme-value Pearson type III Generalized Pareto Generalized logistic Log normal

0.6 0.4

1 2 3 4 5 6 7

0.2

Cool Wet_Lowland

Cool Wet_Mountain

Warm Dry_Lowland

Warm Extremely Wet_Hill

Cool Extremely Wet_Mountain

Cool Wet_Glacial

Cool Wet_Hill

Cool Extremely Wet_Hill

Cool Extremely Wet_Lake

Cool Extremely Wet_Lowland

Cool Dry_Lowland

Cool Dry_Mountain

0.6 0.4

L-kurtosis

0.2 Cool Wet_Lake

0.6 0.4 0.2 Cool Extremely Wet_Glacial

0.6 0.4 0.2

Cool Dry_Hill

Cool Dry_Lake

0.6 0.4 0.2

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

L-skewness Fig. 4. L-kurtosis against L-skewness for each standardised flow duration curve by Source-of-Flow class and stream order (n = 379). Lines represent various probability distribution functions.

20% of the stations were selected randomly. For each cross validation subset, models were fitted to the remaining 80% of the data by incrementally increasing model complexity. At each increment, predictions were made for the 20% of withheld stations for which the predictive performance of the model was evaluated. The complexity of the final model was determined as the number of predictors that produced the maximum average predictive performance. The fitted model was then used to calculate the catchment characteristic dissimilarities between new sites and all possible substitute gauging stations and the gauging station whose dissimilarity was least was used as the substitute FDC. 4.4.3. Distance matrix A Manhatten distance matrix representing differences between all pairs of FDC linear moments (l2, lca and lkur) was calculated. For

each FDC, a substitute FDC was chosen by finding the paired FDCs with the smallest dissimilarity. These estimated FDCs represent those that would have been estimated for ungaunged sites assuming that the catchment characteristics could explain all differences between all pairs of linear moments. In other words: the best estimation of the FDC that could have possibly been achieved using dissimilarity modelling. This method is referred to as the FDC Distance Matrix method in further analysis.

4.4.4. Random substitution For comparison with the above statistical approaches, a further method for estimating the FDC at an ungauged site was defined. For this method, the unknown FDC at an ungauged site was substituted with a randomly chosen observed FDC. This method was

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Probability distribution function

86

Lake (n = 18)

Lowland (n = 155)

Mountain (n = 43)

ALL (n = 379)

Glacial (n = 10)

Hill (n = 153)

10-2 10-1 100 101 102 103

10-2 10-1 100 101 102 103

P3 NORM LP3 LN GUMBEL GEV EV2 P3 NORM LP3 LN GUMBEL GEV EV2

10

-2

-1

10

0

10

1

10

2

10

3

10

ADC Lake (n = 18)

Lowland (n = 155)

Mountain (n = 43)

ALL (n = 379)

Glacial (n = 10)

Hill (n = 153)

Probability distribution function

P3 NORM LP3 LN GUMBEL GEV EV2 P3 NORM LP3 LN GUMBEL GEV EV2 -2000

0

2000 4000 6000

-2000

0

2000 4000 6000

-2000

0

2000 4000 6000

AIC Fig. 5. Anderson–Darling criteria (ADC) and Akaike Information Criteria (AIC) for different probability distributions by hydrological Source-of-Flow.

used to provide a measure of estimation performance through comparison with the other methods. 4.5. Model testing A jack-knife cross-validation procedure (Efron, 1982) was used to provide a test of each method for estimating the FDC at ungauged sites. For each method, this cross-validation procedure was applied by leaving out all data associated with each of the 379 gauging stations and then estimating the FDC for the left-out gauging station using data from all remaining gauging stations. The results from this procedure produced estimates of each FDC for each method as if that gauging station were an ungauged site (Ganora et al., 2009). These jack-knifed comparisons allowed an assessment of both the robustness and reliability of each method for FDC estimation at ungauged sites (Castellarin et al., 2004). After having calculated these jack-knifed estimated FDCs, three methods were employed for assessing different aspects of correspondence between observed and estimated FDCs. The first method was to calculate root-mean-square-deviance (RMSD) for each FDC estimation method for each FDC.

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u nj u X ðlog Q obs ij  log Q est ij Þ2 RMSDj ¼ t nj i¼1

ð16Þ

RMSDj represents a measure of the overall difference between observed and estimated FDC over all exceedance percentiles for each, j, FDC. Zero values for either observed or estimated flows were excluded from calculation of RMSD. Therefore, nj is the number of percentiles with positive values for both observed and estimated FDCs. For the second method, error bands around each flow percentile were calculated by subtracting log estimated standardised flow from log observed standardised flow for each FDC estimation method for each FDC for each flow percentile. Again, zero values for either observed or estimated flows could not be included in this analysis. These error bands were plotted around the mean FDC calculated over all FDCs. As zero flows could not be included in either the first or second methods for model testing, the third method was designed to assess correspondence between observed and estimated proportion of time for which flows were zero. Observed and estimated proportion of time for which flows were zero was compared. Many of the FDC estimation methods did not permit zero flows to be estimated. We set a threshold of standardised flow at a very low value of 0.0001, below which all estimated standardised flows were set to zero. Only 9% of sites had any standardised observed flows that were below this value. The proportion of sites with zero flows was also 9%. Where possible FDCs were also calculated using parameters fitted directly from the observed FDC data (i.e. with no general-

87

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Logl2 Stepwise

2

3

4

4 0

2

3

4

0

1.5

Linear a1 Stepwise

0.0

0.5

1.0

1.5

1.0

1.5

0.5

1.0

1.5

l2 PRG 0.0 0.2 0.4 0.6 0.8 1.0

0.0 0.2 0.4 0.6 0.8 1.0

Comb3 a1 Stepwise 3 2 1 0

0

1

2

0.0 0.2 0.4 0.6 0.8 1.0

0.00.51.01.52.02.53.0

0.0 0.2 0.4 0.6 0.8 1.0

4

Linear a1 PRG

0.0

l2 Stepwise

Comb3 a1 RandomForest

0.00.51.01.52.02.53.0

3

1.0 0.5

0.0 0.2 0.4 0.6 0.8 1.0

l2 RandomForest

0.0 0.5 1.0 1.5 2.0 2.5 3.0

2

0.5 0.0

0.0 0.2 0.4 0.6 0.8 1.0

1

0.0

0.0

0.0

0.5

0.5

1.0

1.0

1.5

Linear a1 RandomForest

1

1.5

1

0

0

1

1

2

2

3

3

4 3 2 1 0 0

Predicted

Logl2 PRG

4

Logl2 RandomForest

3

Comb3 a1 PRG

0.0 0.5 1.0 1.5 2.0 2.5 3.0

Fitted Fig. 6. Fitted versus jack-knifed predictions of various example parameters used to describe the shape of FDCs calculated using different methods. PRG refers to the parameterise then regionalise then generalise method.

isation or regionalisation). These fitted FDCs represent those derived using linear models and probability distribution functions that would have been achieved assuming a perfect method for generalising all required parameters. In other words: the best estimation of the FDC that could have possibly been achieved using PG methods. These FDC estimates are referred to as Fitted in further analysis. 5. Results 5.1. Parameterise then generalise 5.1.1. Linear regression models Summary statistics describing the fitted performance of linear model formulations (Eqs. (1)–(6)) to each FDC were compared (Fig. 3). Adjusted r2 increased consistently across FDCs as terms were added to model formulations. AIC was reduced as terms were added across the suite of model formulations. Reduced AIC for the Comb3 model in comparison with the Squared model indicated that adding

a cubic term provided a more appropriate formulation of more FDCs, and therefore a better fitting formulation for the FDCs. This finding was supported by a reduction in residual standard error for 94% of the FDCs when modelled using the Comb3 model rather than the Squared model. AIC was reduced most when a stepwise procedure was used because this procedure is designed to minimise AIC. Models resulting from stepwise reduction of Eq. (6) for each observed FDC showed that there was some between-FDC variation in the combination of terms that best fitted the shape of the log-normal FDC (Table 7). The most complicated stepwise models included all five possible terms, however many reduced models excluded the log term. Table 7 showed no systematic patterns in the inclusion of each of the terms with Source-of-Flow categories. These results indicated that most, but not all, of the 379 observed FDCs were better fitted by more complicated linear models. The greatest reduction in AIC was made when squared and cubic, but not log, terms were used together. The Comb3 model (Eq. (5)) was therefore identified as a candidate national-scale model from amongst the available linear models (Eqs. (1)–(5)).

88

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Table 8 Summary statistics for mixed-effects models of log standardised flow against normal-reduced exceedance percentile. Each model has a different structure for the random effects. Values in brackets indicate standard errors. Parameter

Symbol

Value Model 0

Model 1

Model 2

Model 3

Model 4

Model 5

0.178 (0.0310) 0.386 (0.0004)

0.141 (0.0005) 0.352 (0.0297)

Fixed effects Intercept Slope

a0 a1

0.143 (0.0082) 0.370 (0.0288)

0.135 (0.0130) 0.354 (0.0278)

0.144 (0.0071) 0.370 (0.0282)

0.143 (0.0073) 0.353 (0.0277)

Random effects Source-of-Flow Intercept Slope Corr

s0 s1 s

0.025 0.109 0.874

0.026 0.118 0.982

0.030 0.106 0.634

0.031 0.118 0.495

t0 t1 t e

0.046 0.094 0.136 0.216 82805.866 82708.331 41411.933

0.099

Stream order Intercept Slope Corr Residual AIC BIC Log likelihood

0.131 0.100

0.225 52929.033 52853.173 26471.517

0.218 79265.201 79189.340 39639.600

Table 9 Rank of predictor variable importance for each random forest.

LogusAnRainVar usAvTCold usAvTWarm LogusRainDays10 LogusRainDays50 usPET LogusCatElev usLake LogusCalc usHard usParticleSize usAveSlope usSolarRadWin usSolarRadSum LogCATCHAREA OrderNumeric

0.237 16232.443 16167.420 8122.222

0.252 29563.896 29607.245 14777.948

0.238 13367.834 13324.485 6687.917

Nearest Gauge REC

Random Site

Dissimilarity

Distance Matrix

10

Response parameter

9

Linear

Comb3

Q

Log Q

a0

a1

a0

a1

a2

a3

l2

lca

l1

l2

lca

1 7 5 3 12 2 4 7 13 15 7 10 5 11 14 16

14 5 4 8 12 7 1 10 13 11 2 9 6 3 15 16

2 12 5 4 8 7 6 16 11 14 15 9 1 3 10 13

13 2 4 11 10 12 9 7 16 7 1 3 5 6 15 14

12 4 9 14 8 15 4 11 10 3 7 2 6 1 16 13

1 8 5 3 11 3 5 14 7 10 12 13 9 2 16 15

10 6 5 4 14 8 1 12 16 12 2 3 7 9 11 15

3 8 11 5 12 9 1 4 15 14 2 6 13 9 7 16

11 4 8 12 7 6 1 3 13 14 2 10 5 8 16 15

13 4 11 7 6 9 1 3 14 12 2 9 8 5 15 16

13 4 6 2 8 10 6 15 8 12 3 11 5 1 14 16

8

Log (Catchment area (m2))

Predictor variable

0.126

7 6

10 9 8 7 6

Different class Same class

6

5.1.2. Probability distribution functions There were strong interrelations between the L-moments calculated from each FDC (Fig. 4). The appropriateness of various probability functions can be assessed by comparing empirical patterns between L-moments and theoretical relationships defined by the probability distribution functions (Hosking and Wallis, 1997). This is because probability distribution functions are defined by relationships between L-moments (Eqs. (7)–(10)) or log transformations of these data (Eqs. (11)–(13)) (Table 5: Laio et al., 2009). For example, a relatively small change in L-kurtosis with changes in L-skewness are indicative of a Pearson type III (P3: Eq. (10)) distribution. A greater rate of change or more curvature, as in Fig. 4, are indicative of a generalised extreme value (GEV: Eq. (9)) or a log Pearson Type III (LP3) distribution respectively (Ganora et al., 2009). Visual inspection of Fig. 4 indicates that a relatively tight relationship between L-kurtosis and L-skewness for the 379 observed FDCs, with different river types (as defined by Source-ofFlow) occupy different locations along this relationship. There was a lack of systematic differences in the relationship between L-kurtosis and L-skewness across stream order. Both AIC and ADC values indicated that the Gumbel (Eq. (9)) and normal (NORM: Eq. (10)) distributions were the least appropri-

7

8

9

10

6

7

8

9

10

Log (Catchment area (m2)) Fig. 7. Comparison of catchment areas, and source of flow category for FDCs being substituted (x-axis) and substitute FDCs (y-axis) chosen using various methods.

ate distributions for the 379 observed FDCs (Fig. 5). AIC values indicated that FDCs from rivers predominately fed by lakes were equally well described by all distributions other than the NORM distribution. However, ADC distributions indicated that these FDCs were better described by the LP3 and GEV distributions. When all sites were considered, regardless of Source-of-Flow class, there was great overlap in AIC values calculated for the Frechet (EV2), GEV, LN, LP3 and P3 distributions. There was less overlap in the distribution of ADC, with the LP3 and GEV exhibiting the lowest spread of ADC values. The LP3 and GEV distributions reduced ADC more than any other distribution for 49% and 25% of the observed FDCs respectively. Further analysis therefore concentrated on the LP3 and GEV probability distribution functions, as both AIC and ADC values indicated that these were candidates for probability functions from which generalised FDCs could be generated (Fig. 5). PG and PRG of linear moments calculated from raw standardised FDC data was implemented to enable calculation of GEV

89

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Random Gauge Dissimilarity Distance matrix Nearest REC Mixed Fitted Mixed LP3 Reg Then Gen

Model description

LP3 Step LP3 RandomForest Fitted LP3 GEV Reg Then Gen GEV Step GEV RandomForest Fitted GEV Comb3 Reg Then Gen Comb3 Step Comb3 RandomForest Fitted Comb3 Linear Reg Then Gen Linear Step Linear RandomForest Fitted Linear 0.0

0.2

0.4

0.6

0.8

1.0

1.2

Root mean square deviance in predicted log standardised flow Fig. 8. Root-mean-square-deviance in predicted log standardised flow for all sites (n = 379) across all flow percentiles (n = 1001 unless zero flow are either predicted or observed) using various estimation methods. Reg Then Gen refers to parameterise then regionalise then generalise strategies, where Linear, Comb3, GEV or LP3 refers to the parameterisation method.

Warm Dry

Warm Wet

Warm Extremely Wet

Cool Dry

Cool Wet

Cool Extremely Wet

7 6 5 4 3

Stream Order

2 1

7 6 5 4 3 2 1 0.0

0.5

1.0

0.0

0.5

1.0

0.0

0.5

1.0

Root mean square deviance in predicted log standardised flow Fig. 9. Root-mean-square-deviance in predicted log standardised flow across all flow percentiles (n = 1001 unless zero flow are either predicted or observed) by climate category and stream order using the GEV Random Forest method.

distributions, and linear moments calculated from the log transformed standardised FDC data was implemented to enable calculation of LP3 distributions.

5.1.3. Generalisation Jack-knifing procedures were used to re-calculate each parameter required for: (a) the Linear model; (b) the Comb3

90

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

model; (c) L-moments of raw observed FDC data; and (d) L-moments of logged observed FDC data. Each parameter was calculated from catchment characteristics (Table 6) using two jack-knifed regression methods: (a) stepwise linear regression with two-way interactions; and (b) random forests. Comparison between fitted and jack-knifed calculated parameter values showed that parameter values for more sites were more accurately calculated using random forests, but that extreme values were more accurately calculated by linear stepwise models. This is because random forests are unable to predict outside of the range of observed values (Fig. 6). Variable importance was calculated for each explanatory variable for each response that was generalised using random forests.

Results indicated that there were differences in importance between explanatory variables within each response (Table 9). There were some consistent patterns between modelled responses, for example, upstream catchment elevation (LogusCatElev) had the highest importance for five of the eleven response variables and particle size in the upstream catchment (usParticleSize) had the second highest importance for six of the eleven responses. However, variation in rainfall (LogusAnRainVar), summer radiation (usSolarRadSum) and potential evapotranspiration (usPET) also had high importance for some responses. This indicates that different FDC characteristics may be linked with different catchment characteristics and that there may be considerable interactions between explanatory predictors within these random forests models.

Fitted Linear

Linear RandomForest

Linear Step

Linear Reg Then Gen

Fitted Comb3

Comb3 RandomForest

Comb3 Step

Comb3 Reg Then Gen

Fitted GEV

GEV RandomForest

GEV Step

GEV Reg Then Gen

Fitted LP3

LP3 RandomForest

LP3 Step

LP3 Reg Then Gen

Fitted Mixed

Mixed

Nearest REC

Distance matrix

Dissimilarity

Random Gauge

1 0 -1 -2

1 0 -1

Observed Log (Q/Qbar) - predicted Log (Q/Qbar)

-2

1 0 -1 -2

1 0 -1 -2

1 0 -1 -2

80th percentile 90th percentile 95th percentile Mean FDC

1 0 -1 -2

0

20

40

60

80 100 0

20

40

60

80 100

Percent of time that flow is not exceeded Fig. 10. Observed minus predicted standardised flow (n = 379) at each flow percentile (n = 1001) plotted around the averaged observed FDC. Note that points with either zero observed or predicted flows cannot be plotted. Reg Then Gen refers to parameterise then regionalise then generalise strategies, where Linear, Comb3, GEV or LP3 refers to the parameterisation method.

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

5.2. Parameterise then regionalise then generalise Comparison between fitted parameter values and those calculated from the jack-knifed estimations after regionalisation, showed that generalisation by averaging within Source-of-Flow groups did not result in improved parameter predictions (Fig. 6). This suggests that, although separating the analysis by river type may have isolated FDCs into groups with more consistent relationships between FDC shape parameters and catchment characteristics, there was still considerable within group variation in parameter values (Fig. 6). 5.3. Parameterise and generalise together The appropriateness of several mixed-effects models for estimating logged standardised Q as a function Z was assessed by comparing model summaries. Six mixed-effects models, each with a different formulation of the random-effects, are shown here as examples (Table 8). AIC, BIC (Bayesian Information Criterion: Schwarz, 1978) and model residuals indicated that inclusion of stream order nested within Source-of-Flow for both terms in Eq. (1) created improvements in explanatory power in comparison with employing either of these grouping factors in isolation. Likelihood ratio tests showed that Model 0 (Eq. (14)) had significantly greater explanatory power than any of the reduced models (p < 0.001). Comparisons between the fixed and random effects for Model 0 indicated that there were considerable variations in a0 and a1 between stream orders within Source-of-Flow groups (Table 8). These results also indicated differences in model performance with changes in the formulation of random-effects, with Model 0 outperforming other formulations. P-values for Wald tests on individual parameters (intercept: F = 1106, p < 0.0001, slope: F = 164, p < 0.001) showed that both fixed effects were significant. Model 0 was therefore used as an example PGT strategy in further analysis. 5.4. FDC substitution Fig. 7 indicates the degree to which each FDC-substitution method was able to select substitute FDCs from gauging stations with matching catchment characteristics in terms of Source-ofFlow category (Table 1) and catchment area. The Nearest in REC class method identified substitute FDCs by selecting the FDC from the nearest gauging station from within the same Source-of-Flow category. Since no additional criteria were used some FDCs were substituted with FDCs whose gauging stations had very different catchment areas. The ‘‘Distance matrix method’’ identified pairs of FDCs with the most similar L-moments. Spread away from the one-to-one line and the prevalence of FDCs paired with FDCs from gauging stations with different Source-of-Flow classes indicates that similar FDC shapes can result from different catchment characteristics. The Dissimilarity method produced a tight fit around the one-to-one line and many FDCs were paired with FDCs from gauging stations with the same Source-of-Flow classes. This was because the dissimilarity method is designed to identify differences in FDCs that can be explained by differences in catchment characteristics. 5.5. Model testing RMSD between observed and estimated log standardised Q was calculated from various FDC estimation methods (Fig. 8). Results from fitted models indicated that, the Comb3 model and LP3 distribution were better able to represent the observed FDCs in comparison with both the Linear model and the GEV distribution. Results suggested that, given perfect generalisation, the LP3 model would

91

have best matched the observed FDCs in comparison with the other methods used to apply a PG strategy. However, after generalisation, RMSD was lower for the Linear model and GEV distribution than for the Comb3 model and LP3 distribution respectively, regardless of generalisation method. This indicated that, although the Comb3 model and LP3 distribution were better able to describe the shapes of the observed FDCs, generalisation of the parameters required for these methods was less reliable than was the case for both the Linear model and the GEV distribution. For the PG strategies, RMSD was reduced most when parameters were generalised using random forests, rather than stepwise linear models, regardless of FDC parameterisation method. Generalisation using random forests also out-performed the PRG strategy, regardless of FDC parameterisation method. RMSD was slightly more for the mixed-effects model than for the PG strategies that were generalised using random forests. However, the mixed-effects model required less input information (only Source-of-Flow and stream order) as opposed to the random forest methods which took all variables listed in Table 6 as input. RMSD for the mixed-effects method was greater than that for the dissimilarity method and the nearest in REC class method. RMSD was not the same across river sizes or catchment types (e.g., Fig. 9). For the GEV Random Forest method, RMSD was generally less for larger rivers and rivers in wetter climates. Although inequality in sample size makes it difficult to draw firm conclusions, these results indicate that links between flow variability and catchment characteristics are less consistent for smaller rivers and for catchments with drier climates. Errors bands were calculated across percentiles and plotted in relation to the average of the 379 observed FDCs. Results indicated that error was not the same across exceedance percentiles (Fig. 10). Estimated standardised FDCs exhibited wider error bands at lower flows than at higher flows when plotted in log space. Results agreed with calculated RMSD in that there was less error for the fitted Comb3 model in comparison with either the fitted LP3 or fitted GEV models. This was even the case for very low flows, where the fitted Comb3 model still performed well. The GEV model tended to overestimate flows in the medium range and underestimate lower flows, whereas the LP3 method did not exhibit systematic over- or under-estimation. Errors increased considerably after jack-knifed generalisation for all four methods that employed PG strategies. This was particularly the case when linear stepwise regression with two-way interactions was used to generalise the Comb3 model, which exhibited nearly as much error as choosing to substitute with a random gauge. Results suggested that, in this case, random forests produced smaller errors in comparison to stepwise linear regression. Errors for the methods used to employ a PRG strategy were an improvement over replacement with an FDC from a random gauge, but errors were still considerable, especially at lower flows. This shows that averaging within Source-of-Flow classes only captures limited between-site variation in the shapes of FDCs. This finding matches well with the comparison of various formulations of the random-effects within mixed-effects models, which showed considerable improvements when stream order was nested within Source-of-Flow (Table 8). After jack-knifing, errors associated with the mixed-effects predictions were comparable with the other methods (Fig. 10). However, the mixed-effects model tended to overestimate low flows in a few FDCs with steep slopes. This indicates both a departure from a linear relationship in log-normal space, and that information on just Source-of-Flow together with stream order, could be used to explicate between-FDC variation for many, but not all sites. The FDC distance matrix method represents the best results that would be achieved if a model explaining all differences in

92

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Fitted Linear

Linear RandomForest

Linear Step

Linear Reg Then Gen

Counts

60 40

361 20 0 Fitted Comb3

Comb3 RandomForest

Comb3 Step

Comb3 Reg Then Gen

244

60 40 123 20

Estimated percentage of time with zero flow

0 Fitted GEV

GEV RandomForest

GEV Step

GEV Reg Then Gen

62

60 40 31 20 0 Fitted LP3

LP3 RandomForest

LP3 Step

LP3 Reg Then Gen

16

60 40 9 20 0 Fitted Mixed

Mixed

Nearest REC

Distance matrix

5

60 40 3 20 0 Dissimilarity

Random Gauge

2

60 40 1 20 0 0

20

40

60

0

20

40

60

Observed percentage of time with zero flow Fig. 11. Observed versus predicted (using various methods) proportion of zero flows for 379 FDCs.

FDC L-moments with differences in catchment characteristics was employed within a dissimilarity model (i.e. a perfect dissimilarity model). This is equivalent to the fitted performance for the dissimilarity model. The difference in errors between FDCs estimated using the distance matrix method and the dissimilarity method can be attributed to the degree to which the dissimilarity model cannot explain differences in FDC L-moments using differences in catchment characteristics. Results indicate that, had the dissimilarity model been able to identify all pairs of stations whose FDCs were most similar using catchment characteristics, this method would have out-performed all other methods for the majority (e.g., 80%) of sites. However, the performance of the particular dissimilarity method implemented here was similar to the relatively simple Nearest in REC method (Fig. 10).

Neither comparison of RMSD nor plotting of errors across percentiles allowed assessment of the ability of each estimation method to estimate zero flows. Where estimated zero flows were those calculated to be below 0.0001. This was because these assessments were made in log space, and zero flows cannot be logged. The number of observed zero flow days for each FDC was compared with that calculated by the various methods (Fig. 11). There was strong correspondence between observed zero flow days and those calculated using the Fitted Comb3 model, Fitted LP3 distribution, and to a lesser extent the Fitted GEV distribution. A considerable loss in performance was evident after jack-knifed generalisation, with little difference in the performance of the random forest method in comparison to the stepwise linear method. However, the generalised GEV and LP3 methods as well as the Nearest in REC and

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

dissimilarity methods were all considerable improvements on previously published results for estimating number of zero flow days in New Zealand (Pearson, 1995). Methods that employed PRG strategies consistently under-estimated the number of zero flow days. This is further evidence to support the presence of within Source-of-Flow class variation in FDC characteristics. All methods that assumed linearity of the FDC in log-normal space, including the mixed-effects model, showed an inability to estimate any zero flows.

6. Discussion This work indicated that, after jack-knifing, the shapes of some FDCs were hard to estimate regardless of methods used. This may be caused by a combination of several factors. First, an attempt was made to use data from catchments that were reasonably natural, but despite these efforts, the dataset may have contained sites whose hydrology regimes were influenced by human activities (e.g., abstraction or water storage). Second, data on catchment characteristics were extracted from available national-scale databases. Improved spatial resolution may have produced more accurate representations of the true catchment characteristics, and therefore produced more accurate estimates of FDCs. Third, some hydrological processes may be difficult to generalise empirically given available information on catchment characteristics, the number of sites, and range of hydrological conditions represented by the dataset. Fourth, it was assumed that the sample of 379 gauging stations was broadly representative of the range of catchments found throughout New Zealand. It was assumed that these gauging stations represented various locations along a continuum of catchments and hydrological regimes found across New Zealand. Therefore a jack-knifing procedure was used to quantify errors for each FDC as if it was an ungauged site, rather than repeatedly holding out sets of sites, as would be done for a cross-validation procedure (Picard and Cook, 1984). This may not have been the case. For example, larger catchments may have been over-represented. Fifth, records of various lengths and covering various time periods were used. This is not ideal since extreme events are more likely to appear in longer records. In this study, consistency of record length and period of record coverage were sacrificed in order to include more FDCs with greater spatial coverage. It was assumed that any errors caused by variation in record length did not bias comparison between FDC parameterisation methods, or generalisation methods. Throughout the analysis we concentrated on the log-normal standardised FDCs. All estimated standardised FDCs could be transformed into units of m3 s1 by raising to the power 10 and multiplying by mean flow. Estimates of mean flow at ungauged sites in New Zealand are available (Woods et al., 2006). Thus our strategies assume that an accurate method of estimating mean flow at ungauged sites was available. However, any errors in estimating mean flow would affect the estimated FDC. Although purely empirical approaches that did not include any attempt to include physical processes were used, many of the patterns that were found were physically meaningful. For example, l2 (the second linear moment) and a1 (the slope of the flow against percentile) both decreased with increasing catchment area and elevation. This is indicative of the less flashy flow regimes expected in larger and lower elevation catchments. Furthermore, of the available explanatory variables, catchment elevation, catchment area, potential evapotranspiration and the number of rain days with rainfall over 10 mm generally had high importance for the random forest models of each parameter. This is consistent with previous broad-scale descriptions of hydrology across New Zealand (e.g., Toebes and Palmer, 1969; Woods et al., 2006). Patterns in flood

93

flows (Pearson and McKerchar, 1989), low flows (Pearson, 1995) and flow variability (Jowett and Duncan, 1990) across New Zealand have been previously linked to precipitation, potential evaporation and catchment area, with more variable flow regimes present in higher elevation, wetter and smaller catchments. Stepwise linear regression for generalisation performed poorly in comparison with random forests. This may be because complex non-linear relationships and high-order interactions exist between FDC parameters and combinations of catchment characteristics. After jack-knifed generalisation, the Comb3 method performed poorly in comparison with the GEV and LP3 methods. This may have been because each parameter was generalised independently, even though they were not independent of each other. In contrast, the mixed-effects model was able to estimate sets of parameters together and therefore provide estimated FDCs that better matched the observed data. For this reason, the performance of the mixed-effects model was similar to that for the other methods despite the mixed-effects model requiring less input information. In this case, stream order nested within Source-of-Flow were used to specify random-effects within mixed-effects models to illustrate that mixed-effects could be used as a PRT strategy. Various alternative grouping factors could have been used to specify the randomeffects. For example, a geology categorisation of the gauging stations could have been added. The REC database also contains categorical data describing the dominant valley landform and landcover (Snelder and Biggs, 2002). Alternative combinations of random-effects were not trailed to avoid further increasing the degrees of freedom within the mixed-effects model, and to avoid debate regarding the nested nature or these categorical variables. 7. Conclusion There are several different strategies and many methods that can be used to estimate FDCs at ungauged sites. For parameterise then regionalise strategies, it was found that the combination of parameterisation method and generalisation methods together, rather than either in isolation, was important in determining the overall performance. Results indicated that predictive performance varied between methods and across exceedence percentiles. The mixed-effects approach provided the most parsimonious method for estimating FDC at ungauged sites. A method using the generalised extreme value probability distribution function generalised using random forests was the most accurate method of estimating flow duration curves at ungauged sites across New Zealand. Acknowledgements This research was funded by the New Zealand Ministry of Science and Innovation, Environmental Flows Programme (C01X1004). We thank Maurice Duncan, Eric Sauquet and an anonymous referee for comments on earlier drafts of this manuscript. References Akaike, H., 1973. Information theory as an extension of the maximum likelihood principle. In: Petrov, B.N., Csaki, F. (Eds.), Second International Symposium on Information Theory, Akademiai Kiado, Budapest, pp. 267–281. Arnell, N.W., 1995. Grid mapping of river discharge. J. Hydrol. 167, 39–56. Bárdossy, A., 2007. Calibration of hydrological model parameters for ungauged catchments. Hydrol. Earth Syst. Sci. 11, 703–710. Booker, D.J., 2010. Predicting width in any river at any discharge. Earth Surf. Proc. Land. 35, 828–841. Booker, D.J., Dunbar, M.J., 2004. Application of Physical HAbitat SIMulation (PHABSIM) modelling to modified urban river channels. River Res. Appl. 20, 167–183. Booker, D.J., Dunbar, M.J., 2008. Predicting river width, depth and velocity at ungauged sites in England and Wales using multilevel models. Hydrol. Process. 22, 4049–4057. Breiman, L., 2001. Random forests. Machine Learning 45, 15–32.

94

D.J. Booker, T.H. Snelder / Journal of Hydrology 434–435 (2012) 78–94

Burn, D.H., 1997. Catchment similarity for regional flood frequency analysis using seasonality measures. J. Hydrol. 202, 212–230. Buytaert, W., Bevan, K., 2009. Regionalization as a learning process. Water Resour. Res. 45, W11419. doi:10.1029/2008WR007359. Castellarin, A., Galeati, G., Brandimarte, L., Montanari, A., Brath, A., 2004. Regional flow-duration curves: reliability for ungauged basins. Adv. Water Resour. 27, 953–965. Chambers, J.M., 1992. Linear models. In: Chambers, J.M., Hastie, T.J. (Eds.), Statistical Models in S. Wadsworth & Brooks/Cole. Cigizoglu, H.K., Bayazit, M., 2000. A generalized seasonal model for flow duration curves. Hydrol. Process. 14, 1053–1067. Clapcott, J.E., Young, R.G., Goodwin, E.O., Leathwick, J.R., 2010. Exploring the response of functional indicators of stream health to land-use gradients. Freshw. Biol. 55, 2181–2199. Clapcott, J.E., Collier, K.J., Death, R.G., Goodwin, E.O., Harding, J.S., Kelly, D., Leathwick, J.R., Young, R.G., 2011. Quantifying relationships between land-use gradients and structural and functional indicators of stream ecological integrity. Freshw. Biol.. doi:10.1111/j.1365-2427.2011.02696.x. Clausen, B., Young, A.R., Gustard, A., 1994. Modelling the impact of groundwater abstractions on low-river flow. In: Seuna, P., Gustard, A., Arnell, N.W., Cole, G.A. (Eds.), FRIEND: Flow Regimes from International Experimental and Network Data, IAHS Publication No. 221, pp. 77–85. Crawley, M.J., 2002. Statistical Computing: An Introduction to Data Analysis Using S-Plus. Wiley, Chichester, UK, 761pp. Croker, K.M., Young, A.R., Zaidman, M.D., Rees, H.G., 2003. Flow duration curve estimation in ephemeral catchments in Portugal. Hydrol. Sci. J. 48, 427–439. Cutler, D.R., Edwards, T.C., Beard, K.H., Cutler, A., Hess, K.T., Gibson, J., Lawler, J.J., 2007. Random forests for classification in ecology. Ecology 88, 2783–2792. Efron, B., 1982. The Jackknife, The Bootstrap and Other Resampling Plans. Society for Industrial and Applied Mathematics, Philadelphia, PA. Fennessey, N.M., Vogel, R.M., 1990. Regional flow-duration curves for ungauged sites in Massachusetts. J. Water Resour. Plann. Manage. (ASCE) 116, 531–549. Ferrier, S., Manion, G., Elith, J., Richardson, K., 2007. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Divers. Distrib. 13, 252–264. Franchini, M., Suppo, M., 1996. Regional analysis of flow duration curves for a limestone region. Water Resour. Manage. 10, 199–218. Ganora, D., Claps, P., Laio, F., Viglione, A., 2009. An approach to estimate nonparametric flow duration curves in ungauged basins. Water Resour. Res. 45, W10418. doi:10.1029/2008WR007472. Gelman, A.G., Hill, J., 2007. Data Analysis Using Regression and Multilevel/ Hierarchical Models. Cambridge University Press, Cambridge, 625pp. Holmes, M.G.R., Young, A.R., Gustard, A., Grew, R., 2002. A region of influence approach to predicting flow duration curves within ungauged catchments. Hydrol. Earth Syst. Sci. 6, 721–731. Hosking, J.R.M., 1990. L-moments: analysis and estimation of distributions using linear combinations of order statistics. J. Roy. Statist. Soc. B (Methodol.) 52, 105–124. Hosking, J.R.M., Wallis, J.R., 1997. Regional Frequency Analysis. Cambridge University Press, Cambridge, UK. Isik, S., Singh, V.P., 2008. Hydrologic regionalization of watersheds in Turkey. J. Hydrol. Eng. 13, 824–834. Jowett, I.G., Duncan, M.J., 1990. Flow variability in New Zealand rivers and its relationship to in-stream habitat and biota. New Zeal. J. Marine Freshwater Res. 24, 305–317. Laio, F., Di Baldassarre, G., Montanari, A., 2009. Model selection techniques for the frequency analysis of hydrological extremes. Water Resour. Res. 45, W07416. doi:10.1029/2007WR006666. Leathwick, J.R., Snelder, T., Chadderton, W.L., Elith, J., Julian, K., Ferrier, S., 2011. Use of generalised dissimilarity modelling to improve the biological discrimination of river and stream classifications. Freshw. Biol. 56, 21–38. LeBoutillier, D.V., Waylen, P.R., 1993. A stochastic model of flow duration curves. Water Resour. Res. 29, 3535–3541. Legendre, P., Lapointe, F.J., Casgrain, P., 1994. Modeling brain evolution from behavior: a permutational regression approach. Evolution 48, 1487–1499. Lichstein, J.W., 2007. Multiple regression on distance matrices: a multivariate spatial analysis tool. Plant Ecol. 188, 117–131. McMahon, T.A., 1993. Hydrologic design for water use. In: Maidment, D.R. (Ed.), Handbook of Hydrology. McGraw-Hill, New York.

Mimikou, M., Kaemaki, S., 1985. Regionalization of flow duration characteristics. J. Hydrol. 82, 77–91. Nathan, R.J., McMahon, T.A., 1992. Estimating low flow characteristics in ungauged catchments. Water Resour. Manage 6, 85–100. Pearson, C.P., 1995. Regional frequency analysis of low flows in New Zealand rivers. J. Hydrol. (NZ) 30, 53–64. Pearson, C.P., McKerchar, A.I., 1989. Flood estimation – a revised procedure. Trans. Inst. Profess. Eng. New Zeal. 16 (2/CE), 59–65. Picard, R., Cook, D., 1984. Cross-validation of regression models. J. Am. Statist. Assoc. 79, 575–583. Pinheiro, J.C., Bates, D.M., 2000. Mixed-Effects Models in S and S-Plus. SpringerVerlag, New York, p. 541. Quimpo, R.G., Alejandrino, A.A., McNally, T.A., 1983. Regionalised flow duration curves for Philippines. J. Water Resour. Plann. Manage. (ASCE) 109, 320–330. Sadri, S., Burn, D.H., 2011. A fuzzy C-means approach for regionalization using a bivariate homogeneity and discordancy approach. J. Hydrol.. doi:10.1016/ j.jhydrol.2011.02.027. Sauquet, E., Catalogne, C., 2011. Comparison of catchment grouping methods for flow duration curve estimation at ungauged sites in France. Hydrol. Earth Syst. Sci. 15, 2421–2435. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Statist. 6, 461–464. Singh, R.D., Mishra, S.K., Chowdhary, H., 2001. Regional flow-duration models for large number of ungauged Himalayan catchments for planning microhydro projects. J. Hydrol. Eng. 6, 310–316. Smakhtin, V.U., 2001. Low flow hydrology: a review. J. Hydrol. 240, 147–186. Smakhtin, V.Y., Hughes, D.A., Creuse-Naudine, E., 1997. Regionalization of daily flow characteristics in part of the Eastern Cape, South Africa. Hydrol. Sci. J. 42, 919– 936. Snelder, T.H., Biggs, B.J.F., 2002. Multi-scale river environment classification for water resources management. J. Am. Water Resour. Assoc. 38, 1225–1240. Snelder, T.H., Hughey, K.F.D., 2005. On the use of an ecological classification to improve water resource planning in New Zealand. Environ. Manage. 36, 741– 756. Snelder, T.H., Cattaneo, F., Suren, A.M., Biggs, B.J.F., 2004. Is the river environment classification an improved landscape-scale classification of rivers? J. North Am. Benthol. Soc. 23, 580–598. Snelder, T.H., Woods, R., Biggs, B.J.F., 2005a. Improved eco-hydrological classification of rivers. River Res. Appl. 21, 609–628. Snelder, T.H., Leathwick, J.R.L., Dey, K.L., 2005b. Definition of the Multivariate Environment River Classification: Freshwater Environments of New Zealand. NIWA Client Report: CHC2005-049. Department of Conservation, New Zealand. Snelder, T.H., Lamouroux, N., Leathwick, J.R., Pella, H., Sauquet, E., Shankar, U., 2009. Predictive mapping of the natural flow regimes of France. J. Hydrol. 373, 57–67. Snijders, T.A.B., Bosker, R.J., 1999. Multilevel Analysis: An Introduction to Basic and Advanced Multilevel Modelling. Sage, London, p. 266. Toebes, C., Palmer, B.R., 1969. Hydrological Regions of New Zealand. New Zealand Ministry of Works, Miscellaneous Hydrological Publication No. 4, p. 45. Vandewiele, G.L., Elias, A., 1995. Monthly water balance of ungauged catchments obtained by geographical regionalization. J. Hydrol. 170, 277–291. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics with S, fouth ed. Springer, New York, p. 495. Vogel, R.M., Fennessey, N.M., 1993. L moment diagrams should replace product moment diagrams. Water Resour. Res. 29, 1745–1752. Vogel, R.M., Fennessey, N.M., 1994. Flow duration curves I: new interpretation and confidence intervals. J. Water Resour. Plann. Manage. (ASCE) 120, 485–504. Vogel, R.M., Fennessey, N.M., 1995. Flow duration curves II: a review of applications in water resources planning. J. Am. Water Resour. Assoc. 31, 1029–1039. Wagener, T., Wheater, H.S., 2006. Parameter estimation and regionalization continuous rainfall-runoff models including uncertainty. J. Hydrol. 320, 132– 154. Warnick, C.C., 1984. Hydropower Engineering. Prentice-Hall, Inc., Englewood Cliffs, New Jersey, pp. 59–73. Wiltshire, E.E., 1986. Regional flood frequency analysis I: homogeneity statistics. Hydrol. Sci. J. 31, 321–333. Woods, R.A., Hendrikx, J., Henderson, R.D., Tait, A.B., 2006. Estimating mean flow of New Zealand rivers. J. Hydrol. (NZ) 45, 95–110. Yu, P.S., Yang, T.C., Wang, Y.C., 2002. Uncertainty analysis of regional flow duration curves. J. Water Resour. Plann. Manage. (ASCE) 128, 424–430.