The effects of spatial discretization and model parameterization on the prediction of extreme runoff characteristics

Journal of Hydrology 392 (2010) 54–69

Contents lists available at ScienceDirect

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

The effects of spatial discretization and model parameterization on the prediction of extreme runoff characteristics Rohini Kumar *, Luis Samaniego, Sabine Attinger UFZ Helmholtz-Centre for Environmental Research, Permoserstrasse-15, 04318 Leipzig, Germany

a r t i c l e

i n f o

Article history: Received 31 August 2009 Received in revised form 7 July 2010 Accepted 31 July 2010

This manuscript was handled by Konstantine P. Georgakakos, Editor-in-Chief, with the assistance of Kieran M. O’Connor, Associate Editor Keywords: High and low flow Parameterization Calibration MPR HRU Hydrologic model mHM

s u m m a r y Water resources management in mesoscale river basins requires, among other things, reliable predictions on extreme runoff characteristics such as magnitude and frequency of floods and droughts. Hydrologic models are increasingly used for these prediction purposes. Outputs of these models, however, are sensitive to various factors like the spatial representation of hydrologic processes, the parameterization method, and the type of estimator used for calibration. This study aimed to investigate the possible effects of these factors on extreme runoff characteristics derived from simulated streamflow. For this purpose, lumped and distributed versions of the conceptual mesoscale hydrologic model (mHM) were implemented in 22 German basins ranging in size from 58 to 4000 km2. The distributed mHM version was, in turn, parameterized with hydrological response units (HRU) and multiscale parameter regionalization (MPR) methods. Free parameters of both model versions were calibrated with three objective functions emphasizing high flows, low flows, and a combination of both. Six extreme runoff characteristics were derived from daily streamflow simulations for winter and summer. Results indicated that the model performance evaluated with both daily streamflow and seasonal runoff characteristics was sensitive to the type of estimator, the spatial discretization, and the parameterization method employed. The lumped version exhibited the highest sensitivity to previous factors and the least performance, whereas the opposite behavior was noticed for the distributed version parameterized with the MPR technique. Furthermore, the efficiency of the model parameterized with MPR were higher than that obtained with the HRU parameterization, in particular, when the model was evaluated in locations not used for calibration. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Reliable estimates and predictions of extreme runoff characteristics like magnitude and frequency of floods and droughts are required for water resources management and engineering. Numerous studies in recent years have investigated the expected changes of runoff characteristics in many river basins due to climate and/or land cover changes (e.g. Middelkoop et al., 2001; Lehner et al., 2006; Kundzewicz et al., 2007). In general, results of these studies point that changes in extreme runoff characteristics may be one of the most significant consequences of climate change. Extreme runoff events pose serious risks to human life and entail substantial socioeconomic and environmental damages. Therefore, improving the predictions of extreme runoff characteristics has become one of the major objectives of contemporary hydrology. In general, two approaches have been pursued for quantifying extreme runoff characteristics in gauged basins. The first approach relies on using statistical methods to find parsimonious and robust * Corresponding author. Tel.: +49 235 1972. E-mail address: [email protected] (R. Kumar). 0022-1694/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2010.07.047

data-driven models that relate runoff characteristics with a set of explanatory variables in a given spatial and temporal domain. Examples of this type of analysis can be found in Pilgrim (1987), Samaniego and Bárdossy (2007), and Ouarda et al. (2008). The second approach derives extreme runoff characteristics using continuous streamflow simulated with a hydrologic model (e.g. Booij, 2005; Lehner et al., 2006; Davison and van der Kamp, 2008). This approach has been frequently used during recent years because hydrologic models are based on conceptualizations or physical descriptions of hydrological processes. This feature enables them to adapt to changing environmental conditions such as climate and/or land cover (Bronstert et al., 2002), to account for antecedent conditions of relevant state variables (e.g. soil water content), and to react to extreme events much better than simple data-driven models. Outputs of these hydrologic models are known to be dependent on various factors such as: (1) the spatial representation of hydrological processes, (2) the method used to parameterize a model, and (3) the procedure used to estimate effective model parameters. However, there is still little known about the sensitivity of extreme runoff characteristics to these factors.

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

Many modeling studies have indicated that hydrological response of a basin is sensitive to the spatial heterogeneity of basin physical characteristics such as topography, soil texture and land cover (e.g. Grayson and Blöschl, 2000; Bronstert et al., 2002), as well as to the spatio-temporal variability of meteorological inputs such as precipitation (e.g. Krajewski et al., 1991; Koren et al., 1999a). Therefore, it is expected that spatially distributed hydrologic models outperform their lumped counterpart which do not account for such variabilities. In practice, however, mixed results have been reported in literature. Reed et al. (2004), for example, concluded that ‘‘. . . lumped model outperformed distributed models in more cases than distributed models outperformed the lumped model. . . ”, based on the findings of a comprehensive inter-comparison between several distributed models with lumped model simulations. Other studies (e.g. Perrin et al., 2001; Booij, 2005; Das et al., 2008) have also reported no significant improvement in daily streamflow simulations with increase in model complexity and spatial resolution in both input data and modeling scale. Given such contradictory conclusions and knowing that the growing availability of remote sensing data sets and computational capability will continue to foster the development of distributed hydrologic models, it is not yet clear how such models would help to improve streamflow simulations over their lumped counterparts. One of the major obstacles for the implementation of a distributed hydrologic model in mesoscale river basins is its parameterization (Grayson and Blöschl, 2000). This is due to the large number of unknown parameters that need to be estimated through the calibration process. The high dimensionality of the parameter search space not only imposes a severe challenge for the existing optimization algorithms (Pokhrel et al., 2008), but also tends to increase the predictive uncertainty of model outputs (Seibert, 2003) due to the equifinality of feasible solutions (Beven, 2006). To deal with this problem, several parameterization methods have been developed in the recent past. The most commonly used methods rely on reducing the spatial variability of model parameters by grouping grid cells into somehow homogeneous regions called ‘‘Hydrologic Response Units” (HRU) (Leavesley et al., 1983), and thereby reducing the dimensions of the parameter space by assigning the unique set of model parameters to each HRU type. Recent applications of this approach were reported by Beldring et al. (2003) and Das et al. (2008) to parameterize a spatially distributed version of the HBV model. In these cases soil properties and land cover classes were used to identify individual HRUs. More recently, Blöschl et al. (2008) followed a similar approach to parameterize a grid-based flash flood forecasting model. Reasonably good model performances for streamflow simulations were reported in all these studies. Another possibility to parameterize a distributed hydrologic model is based on a regionalization method, which attempts to link model parameters with catchment characteristics through a priori defined functional relationships. It has been shown in past studies that this approach is not only suitable for facilitating the application of a distributed hydrologic model in ungauged locations, but also for reducing the dimensionality of the optimization problem. Hundecha and Bárdossy (2004), for example, employed this approach to parameterize a semi-distributed HBV model and reported reasonably good performance of the simulated daily streamflow at gauged and assumed ungauged locations (i.e. one which is not the part of model calibration). Pokhrel et al. (2008) followed a similar approach to parameterize the HL-DHM model and reported a significant reduction in the number of calibration parameters (from 858 to 33), without loosing the model predictability. A multiscale parameter regionalization method was recently presented by Samaniego et al. (2010b), which explicitly takes the sub-grid variability of catchment characteristics into ac-

55

count. This method not only reduced the dimensionality of parameter searching space but also facilitates the transferability of model parameters to scales and locations other than those used during calibration, while keeping the model performance within reasonable ranges. To our knowledge, no study has yet been reported in literature that compares the effects of the HRU approach and parameter regionalization technique on model performance. It is also recognized that the parameterization and the performance of a hydrologic model depends on the estimator (i.e. objective function) used during calibration (e.g. Gupta et al., 1998; Wagener and McIntyre, 2005; Oudin et al., 2006). Often, the use of different objective functions results in different model simulations, which may fit specific aspects of the hydrograph at the expense of others. As noticed by Oudin et al. (2006), simulated low flows obtained from TOPMODEL or GR4J models were biased when these models were calibrated using the sum of squared errors as an objective function. On the contrary, both model simulations were biased to high flows when the calibration was performed using the logarithmic transformation of streamflow as objective function. Such sensitivity of the model performance might cause practical problems since a given model is often required for various purposes, for instance in flash flood forecasting and reservoir operation, where efficient predictions for both high and low flow regimes are required. The sensitivity of the estimator on extreme runoff characteristics, however, has not been reported in literature. The objective of this study was to explore the effects of the spatial discretization and model parameterization on the daily streamflow simulations and on extreme seasonal runoff characteristics. The mesoscale hydrologic model (mHM) (Samaniego et al., 2010b) selected for this purpose was applied in many southern German basins at two spatially discretized versions: lumped and distributed. The distributed mHM version was parameterized with two methods: hydrological response units and the multiscale parameter regionalization. We also carried out a sensitivity analysis of the model performance to the type of estimators used for the calibration process. For this analysis, three estimators were formulated emphasizing: high flows, low flows and a combination of both. Split sampling and proxy basin tests were conducted to evaluate the model performance.

2. Study area and data The study was conducted in 22 river basins ranging in size from 58 to 4000 km2 (Fig. 1). They are located in the upper Neckar river basin (Germany). This basin is bordered by the north-western edge of the Swabian-Jura on the right bank of the Neckar river and by the Black Forest on its left bank. The elevation ranges from 250 m to 1013 m with a mean elevation of 550 m a.s.l. Slopes are in general mild with 90% of the area having mild slope varying from 0° to 15°. The climate of the region is moist with mild winters according to the Köppen notation. The long term annual mean precipitation is approximately 900 mm, with the rainiest month being June and the driest one being October. The daily mean air temperature is 8.1 °C, ranging between 19 °C and 27 °C with the coldest and warmest months being January and July, respectively. The required meteorological data to drive the mHM model are precipitation, mean air temperature, and potential evapotranspiration. The daily data from 294 rain gauges and 157 weather stations, which are located within and around the study area, were obtained from LUBW1 and DWD2 for the period from 1979 to 2001. The External Drift Krigging (Ahmed and Marsily, 1987) procedure was 1 2

State Institute for Environmental Protection Baden-Württemberg, Germany. German Meteorological Service.

56

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

2220

172

12 11

N

r ve r ri 7 a k ec

9

Germany

1

10 4

8

Stuttgart

21 5

Ju ra

Fore st

16 13

Gauging station

Sw ab ia n

Blac k

3 15 19

River Basin boundary

14

186

Elevation [m] High : 1013

10

0

km 20

Low : 250

Gauge No. 2

1

3

60

4

5 6 78 9

10

11 121314

15

16 17

200

80

100

21

19 20

18

400

600

800

Upstream area (km2)

22

2000 4000

1000

Fig. 1. Location of the upper Neckar river basin and 22 gauging stations employed in this study. Basins are numbered from smallest to largest drainage area, which is shown in the sub-plot (below) of this figure.

subsequently used to interpolate the point scale precipitation and air temperature data at a spatial resolution of 4 km  4 km (level-2; see Section 3.1) for the distributed hydrological modeling. The terrain elevation was used as an external drift for both meteorological variables. The daily total potential evapotranspiration was estimated with the Hargreaves and Samani method (Hargreaves and Samani, 1985) at the same spatial resolution. The geo-spatial data set required to set up the mHM model comprises of: a digital elevation model (30 m  30 m), a digitized soil map (1:200,000), and a digitized geological map (1:600,000). These data were obtained from LUBW. Land cover scenes were obtained from Landsat (30 m  30 m) for the years 1984 and 1993, which were subsequently classified and grouped into three major classes: forest cover, permeable areas and impervious areas. All these data were finally compiled to a common scale of 100 m  100 m (level-0; see Section 3.1). Weekly leaf area index data were acquired from the Moderate Resolution Imaging Spectroradiometer (MODIS)3 for the period from 2001 to 2008. These data were further processed to estimate long term mean monthly LAI for different land cover classes and kept constant during the modeling period. The time series of daily streamflow data at all 22 gauging stations (Fig. 1), for the period from 1979 to 2001, were obtained from LUBW. Readers may refer to Samaniego et al. (2010b), for more details regarding data processing.

3

http://wist.echo.nasa.gov.

3. Model description and parameterization 3.1. The mesoscale hydrologic model (mHM) A conceptual mesoscale hydrologic model (mHM) successfully tested in more than 30 basins in Germany (Samaniego et al., 2010a) was used in this study. The basic formulation of this model for describing various hydrological processes is based on the well known HBV model (Bergström, 1995; Hundecha and Bárdossy, 2004). mHM also includes a number of new features to enhance its predictability (Samaniego et al., 2010b), namely: canopy interception, two-soil layers in the root zone, soil freezing and thawing processes as well as cell to cell routing, amongst others. Two model versions were formulated (lumped and distributed) to study the effect of spatial discretization. The lumped version of mHM represents the study basin as a single and homogeneous unit, and thereby does not account for the spatial variability of meteorological forcings and basin characteristics. The distributed version, on the contrary, uses grid cells as an elementary hydrological unit to simulate the spatio-temporal evolution of hydrological processes. The simulated processes in mHM comprise: canopy interception, snow accumulation and melting, soil moisture, soil freezing and thawing, evapotranspiration, infiltration, surface and subsurface runoff, percolation, base flow and discharge attenuation as well as flow routing processes. The schematic representation of each component accounted for in mHM is shown in Fig. 2. Readers interested in more details on this model may refer to Samaniego et al. (2010b). The model requires 28 parameters for every grid cell to

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

Meteorological forcings Temperature T

Potential evapotranspiration Ep

Precipitation P

E1 Interception

1

F S Snow-pack

2

E3

E k2-I Layer - I

I

Layer - II

Soil Moisture storage

Impervious cover

x k3

Ik Upper storage

E k2

6

q4 q1

4

q2

5

q3

C Lower storage

This large optimization problem constitutes a daunting task even for the most advanced optimization algorithms available (Pokhrel et al., 2008). Moreover, calibrating the model with such a large number of unknown parameters for every grid cell will lead to model over-parameterization. This problem, in turn, may leads to an increase of the predictive uncertainty of model outputs (Beven, 2006), and thus the reliability of model outputs. In this study, two parameterization schemes were implemented within the distributed mHM version to overcome this problem, namely: hydrological response units and parameter regionalization.

R

M

k -1 k

57

K

Stream storage

7r

Qr

mesoscale hydrological model (mHM) Fig. 2. Schematic representation of different components accounted in the distributed mHM version. Where X = state variable, E = actual evapotranspiration, q = component of runoff, S = snow precipitation depth, R = rain precipitation depth, F = throughfall, I = infiltration capacity, C = percolation, K = gain/loss flux in a leaking cell, Qr = net runoff produced at the outlet of a grid cell.

describe the spatial dynamics of hydrological processes. The brief description of these parameters is provided in Table 1. The distributed version of mHM requires three levels of gridded information. The upper level (level-2) contains the meteorological inputs (drivers), the middle level (level-1) is used to describe the hydrological processes mentioned above, and finally, the lower level (level-0) contains the catchment characteristics. In this study, level-2 and level-1 were set to the equal spatial resolution (4 km  4 km), while the level-0 was set to 100 m  100 m. The temporal resolution was set to hourly time steps to fulfill the constrains of the flow routing process. Hourly streamflow produced at each grid node was aggregated to daily values to be able to compare them with gauged data. The lumped mHM version simulates streamflow at a given catchment outlet on daily time steps using areal average meteorological forcings. It should be noted that the cell to cell routing process of the distributed model version was not included in the lumped version of the mHM model. Instead the runoff generation mechanism of the basin was modeled with a triangular unit hydrograph similar to that of HBV (Bergström, 1995). 3.2. Distributed model parameterization Each of the 28 parameters per modeling cell required in the distributed version of mHM has to be estimated through calibration.

3.2.1. Hydrological response units (HRU) According to Flügel (1995), ‘‘hydrological response units (HRUs) are distributed, heterogeneously structured entities having a common climate, land use, and underlying topo-pedo-geological associations controlling their hydrological dynamics”. The crucial assumption made in this method is that the variation of hydrological processes within an HRU must be small compared with the dynamics in a different HRU. There is no agreed rule for the delineation of HRUs, although it is one of the most commonly used methods for the parameterization of distributed models. In this study, we used the k-means clustering algorithm to group the modeling cells (level-1) into 15 HRUs. Catchment characteristics used for the delineation include elevation, mean slope, aspect, soil texture, types of geological formations, and land cover classes. Climatic indicators (e.g. long term annual total precipitation) in the study area are correlated with elevation, hence only the latter was used in the HRU classification. Catchment characteristics for each grid cell at level-1 were aggregated from the gridded data available at level-0. Each HRU was further assigned to one of the three land cover classes (i.e. forest, impervious and permeable cover) based on the dominant land cover class. Three different allocation modes were used to assign model parameters. The first one allocates a unique value to those calibrated model parameters which do not exhibit spatial variability, regardless of the HRU definition of a given cell. For example, the parameter threshold temperature used for the conversion of precipitation into snow or rain. The second mode assigns parameters which depend only on land cover class, for example those parameters related with the snow accumulation and melting processes. In this case, those parameters were firstly calibrated for each land cover class independently. The effective parameter assigned to a given HRU was then estimated as the weighted combination of them. The fraction of each land cover in every HRU was used as weighting factor. Finally, the third mode assigns parameters which are dependent upon several catchment characteristics, for example those related to soil texture and land cover. In this case, these parameters were independently calibrated to each HRU type, and than assigned to all those modeling cells having similar HRU irrespective to their location. A similar procedure has been followed by Das et al. (2008). A detailed description of these allocation modes for all model parameters is shown in Table 1. The total number of unknown parameters (n), to be estimated through calibration, for the HRU method is n = 4  1 + 6  3 + 18  15 = 292, based on the number of parameters and allocation modes shown in the previous table. Consequently, a significant reduction from NP = 314  28 = 8792 to n = 292 was attained with this method. Here, N denotes the total number of modeling cells at level-1 (i.e. 4 km), and P being the number of mHM parameters per modeling cell. 3.2.2. Multiscale parameter regionalization (MPR) The MPR method for obtaining the spatial fields of model parameters used in this study was successfully employed to parameterize mHM (Samaniego et al., 2010b). This is a two step parameterization procedure. In the first step, the model

58

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

Table 1 List of mHM model parameters. Additionally, the predictor variables and the corresponding upscaling operator: Arithmetic mean (A), Harmonic mean (H), Geometric mean (G) and Majority statistic (M), for the MPR method are shown. Predictor variables include: soil textural data (sand and clay content (%) and root zone depth (mm)), land cover (Forest (F), Impermeable (I) and Permeable (P) cover) and geological characteristics. For the HRU method, parameters range are shown which were assigned by three allocation methods: E = assigned equally to all modeling cells, W = weighted according to the fraction of land cover class of a HRU, and S = being the HRU specific parameter which are assigned to each HRU independently. Param.

b1 b2 b3 b4 b5 bk6

Description

Thickness of waterfilm on the canopy surface (–) Threshold temperature for phase transition snow and rain (°C) Degree day factor during rainless days (mm d1 °C) Rate of increase of the degree-day factor per unit of precipitation (mm d1 °C) Max. degree-day factor reached during rainy days (mm d1 °C) Max. soil moisture content of kth root zone layer (mm)

MPR

HRU

Predictor variable

Upscal. oper.

Min.– Max.

Assign. method

Landcovera –b Land coverb –b Land coverb Soil texture,

A – A – A H

0.1–0.4 1.0–1.0 1.0–4.0 0.1–0.8 4.0–8.0 10.0– 100.0, k=1 100.0– 500.0, k=2

W E W E W S

Land coverc

b7 b8

b9 b10 b11 b12 b13 b14 b15 b16

c d e f g h i j k l m

Soil texture, Land coverd Soil texture,

H H

1.0–5.0 60.0– 80.0

S S

above which the soil is practically impermeable (mm) Shape factor of the gamma distribution that statistically estimates the virtual impermeable area due to frozen soil (–) Antecedent Temperature Index (ATI, proxy for soil temperature) threshold below which unfrozen water content reach its minimum value (K) ATI threshold above which no frozen water exist (K)

Land covere –e

Min. value of unfrozen water content estimated as the fraction of total total water content of first root zone layer (–) Weighting multiplier to estimate ATI from air temp (–) Max. ponding retention in impervious areas (mm) Permanent wilting point, estimated as the fraction of maximum soil moisture content (–) Soil moist. limit above which the actual evapotrans. is equated with the potential evapotransp., estimated as the fraction of maximum soil moisture content (–) Fraction of roots in the first root zone layer (–)



2.0–5.0

S

f

Soil texture

H

S

Soil texturef

H

Soil texturef

H

265.0– 271.0 271.0– 274.0 0.2–0.5

– Landcover Soil texture, Land coverg Soil texture, Land coverg

 A

0.1–1.0 0.1–1.0

S W

H

0.05–0.3

S

H

0.8–1.0

S

g

A

S

Max. water holding capacity of the unsat. zone (mm) Fast recession constant (d) Slow recession constant (d)

Soil texture, land coverh Slopeh Soil textureh

H A A

b21 b22

Exponent that quantifies the degree of nonlinearity of the cell response (–) Effective percolation rate (d)

Soil texturei Soil texturej

A A

b23

Baseflow recession rate (d)

Geologyk

M

k

Geology

M

Length, slope and land cover along drainage path within cell Length, slope and land cover of river reachl Slope of river reachl Aspectm

G

1.0–3.0

S



0.1–10.0

E

 A

0.0–0.5 0.9–1.1

E W

b25

Fraction of the groundwater recharge that might be gained or lost either as deep percolation or as intercatchment groundwater flow in nonconservative catchments (–) Duration of triangular unit hydrograph (h)

b26

Muskingum travel time parameter (h)

b27 b28

Muskingum attenuation parameter (–) Aspect correction factor of the PET (–)

Dickinson (1984). Hundecha and Bárdossy (2004) and Götzinger and Bárdossy (2007). Zacharias and Wessolek (2007). Uhlenbrook et al. (2004). Koren et al. (1999b). Patterson and Smith (1981). Kutilek and Nielsen (1994). Booij (2005) and Uhlenbrook et al. (2004). Booij (2005) and Götzinger and Bárdossy (2007). Liang et al. (1994). Hundecha and Bárdossy (2004) and Le Moine et al. (2007). Tewolde and Smithers (2006). Shevenell (1999).

Land cover

S

b18 b19 b20

b24

a

Parameter that determines the relative contribution of rain or snowmelt to runoff (–) Critical value of soil ice content in first root zone layer

0.1–0.5, for F} 0.8–1.0, for I} 0.5–0.8, for P} 1.0–40.0 1.0–10.0 5.0– 100.0 0.0–1.0 10.0– 100.0 10.0– 10000.0 0.75–1.5

b17

b

S

W

S S S S S S S

59

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

These characteristics were derived from daily streamflow simulations at a seasonal time scale (i.e. summer and winter). The peak over threshold method (Stedinger et al., 1993) and the truncation level method (Smakhtin, 2001) were used to derive the runoff characteristics related to high flows and low flows, respectively. The threshold and the truncation levels were set to the 95th percentile (F0.95) and the 10th percentile (F0.10) of the daily observed streamflow, respectively. These levels were kept fixed for both winter and summer throughout the modeling period. For more details regarding the estimation these runoff characteristics refer to Samaniego (2003).

parameters at the finer resolution (level-0) are regionalized through linear or nonlinear transfer functions by linking them with basin characteristics. The linkage is done through global parameters. Process understanding and empirical evidence (e.g. pedotransfer functions) are used to define these functional relationships. Subsequently, in the second step effective parameter value at the modeling scale (level-1), required to describe the hydrological processes accounted in mHM, are linked with their corresponding ones at the finer resolution (level-0) through an upscaling operator. The general form of these two steps can be summarized as follows

bpi ðtÞ ¼ Op hbpj ðtÞ 8j 2 iii

3.4. Parameter identification

ð1Þ

Good sets of model parameters for both mHM versions were identified using Dynamically Dimensioned Search (DDS) algorithm (Tolson and Shoemaker, 2007). This algorithm is a robust and parsimonious global optimization algorithm that mimics the manual calibration by dynamically and probabilistically reducing the number of parameters during the search procedure. DDS is a single objective optimization algorithm and requires the definition of the overall objective function (or estimator). We tested three different variations of the objective function to assess the sensitivity of the model performance. Each objective function was formulated based on the daily streamflow with an emphasis on high flows, low flows, and a combination of both. Additionally, seasonal extreme runoff characteristics were used as a penalty term to improve the model predictions at this temporal scale. The overall objective function U to be minimized is

where

bpj ðtÞ ¼ fp ðuj ðtÞ; cÞ

ð2Þ

Here p = 1, . . . , P with P = (28) denoting the number of distributed model parameters. uj denotes a v-dimensional predictor vector for cell j at level-0 which is contained by cell i at level-1. Ophii denotes the type of operator applied for the regionalization of the parameter p. c is a s-dimensional vector of global parameters to be calibrated. By establishing such relationships, the calibration algorithm finds good solutions for the global parameters (s = 64) instead of estimating parameters for each grid cell independently, subject to the condition that model predictability is not lost. This is advantageous, because a great reduction in complexity can be achieved since NP  s. A summary of the catchment characteristics and upscaling operators used in the MPR method is shown in Table 1. It is worth noting that the number of global parameters (c) needed to be estimated through calibration remain constant regardless of the changes in modeling scales as shown in the study by Samaniego et al. (2010b). Interested readers may also refer to the same study for more detailed description on the MPR method.

U¼

! C Y ð2  Hc Þ ð1  EÞ

ð3Þ

c¼1

where C denotes the total number of penalty functions and Hc is the Pearson correlation coefficient between observed and simulated runoff characteristic c. In particular, the specific volume of high flows (Q1) and the cumulative specific deficit (Q4) were used as penalty for the high and low flow calibrations respectively, whereas both characteristics were used in the combined calibration case. E is the Nash–Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970), which is formulated as

3.3. Seasonal high and low flow characteristics In general, modeling of extreme runoff characteristics constitutes a challenging task because of their highly skewed probability density functions (Samaniego and Bárdossy, 2007). Consequently, to be able to predict them is a strong test for the performance of a hydrologic model. Six high and low flow characteristics, based on the previous work of Samaniego (2003), were used in this study. They are: specific volume (Q1), total duration (Q2) and frequency (Q3) of high flows, as well as cumulative specific deficit (Q4), total drought duration (Q5), and maximum drought intensity (Q6) of low flows. The graphical representation of these characteristics is depicted in Fig. 3.

P 8 b2 > Pt ðqt  qt Þ2 for high flows > ¼ 1  E Q >  > ðq qt Þ > t t < P d E¼ ðlog qt  log q t Þ2 > Elog Q ¼ 1  Pt for low flows > > ðlog qt log qt Þ2 > t > : 1 EQ log Q ¼ 2 ðEQ þ Elog Q Þ for combined high and low flows ð4Þ αQ 2

αQ 3

60

αQ 1

F0.95

45

q(t)

αQ 5 30

15

F0.10

αQ 4

αQ 6

0 0

40

80

120

160

200

Day Fig. 3. Schematic representation of high flow and low flow characteristics: specific volume (Q1), total duration (Q2), and frequency (Q3) of high flows, specific cumulative deficit (Q4), total drought duration (Q5), and maximum drought intensity (Q6). F0.95 and F0.10 are the threshold limits for defining high flow and low flow characteristics, respectively.

60

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

Here, EQ is the NSE between the daily observed (q) and simu^Þ streamflow to emphasize on high flows. Elog Q, on the conlated ðq trary, places emphasis on the matching of recessions and low flow spells. Finally, the third objective function EQlogQ is the equally  and log q are mean values average of the previous two statistics. q of the observed and their logarithmic transformation over the calibration period, respectively. Equifinality is the principle by which an open system can reach a given end state by many potential means. This principle has been recognized as one of the major problems for environmental modeling (Beven, 2006), and constitutes one of the various causes of predictive uncertainty in hydrological modeling. In the context of this study, this principle indicates that an acceptable level of performance (estimated for instance by EQ) can be reached by different sets of model parameters. The predictive uncertainty of model outputs is considerably increased as a result of the equifinality of parameter sets (Seibert, 2003). A Monte-Carlo approach was carried out to assess the variability in model outputs (e.g. streamflow) due to equifinal parameter sets. For this purpose, a total of 200 independent realizations of the DDS algorithm were carried out for each objective function, model version, and parameterization method, respectively. Each DDS realization comprises of at least 200 model simulations, each of them initialized with a different initial solution. The number of model runs per DDS realization was set up proportionally to the number of free calibration parameters. The distributed model with HRU parameterization got the highest number of model runs per DDS realization, whereas the lumped model version got the least. The NSE value of 0.4 for both estimators (i.e. EQ P 0.4and ElogQ P 0.4) was set as a threshold to identify ‘‘good” parameter sets. The above approach to approximate the variability in model outputs due to good sets of model parameters is based on the DDS-approximation of uncertainty method (Tolson and Shoemaker, 2008). This method is more effective and efficient than the simple random search method (Seibert, 2003; Beven, 2006) often reported in literature. The deterministic model prediction for a particular case was estimated as the median of the ensemble of simulations obtained with the corresponding 200 optimized parameter sets. The variability of the model outputs was estimated with the interval between 5% and 95% percentile ranges, hereafter denoted as (P5–P95). It should be noted that this variability bound does not correspond to the 90% confidence interval in a statistical sense.

4. Application Both mHM versions with their respective parameterization methods were set up independently in all investigated river basins (Fig. 1) to simulate the observed streamflow at their respective outlets. Model runs were performed during the period 1980– 2001. This period was split into two parts, 1980–1988 for the calibration phase, and 1989–2001 for the validation one. Six months of data, prior to the calibration period were used as spin-up time to establish reliable initial conditions of the state variables. The simulations during the spin-up time were not accounted in the estimation of the overall objective function during the calibration process. Two tests were conducted in this study to assess the model performance, namely: a differential split-sample test and a proxy basin test. In the split-sample test a part of data were withheld from the model calibration and then were used for the model evaluation. In the proxy basin test, on the contrary, no direct calibration of parameters at a given location was allowed. Here, the model simulations were performed in a given location with parameters transferred from other place at which the model was calibrated. The

proxy basin test was used to assess the robustness of the distributed model parameterization methods. Free parameters of HRU and MPR obtained through calibration in two locations: gauge no. 21 (Horb; area 1120 km2) and 22 (Plochingen; area 4000 km2) were used to simulate the streamflow characteristics in 21 remaining locations. These two basins were used because their drainage area is large enough to cover all delineated HRUs and, hence, provided adequate support for the inference of model parameters. The performance of mHM for the predictions of seasonal extreme runoff characteristics was judged by two statistics: the Root Mean Square Error (RMSE) and the Pearson correlation coefficient (r).

5. Analysis of results and discussion 5.1. Sensitivity of the model efficiency estimated with daily streamflow Nash–Sutcliffe efficiencies based on daily streamflow simulations obtained with 200 optimized parameter sets for each objective function, model version, and parameterization method are depicted in Fig. 4. All simulations shown in this figure were evaluated at the outlet of the study area (gauge no. 22) for both the calibration and validation periods. Trade-offs in model performance of lumped vs. distributed model versions, and of MPR vs. HRU parameterizations for the distributed version for three types of objective functions can be clearly distinguished in this figure. A consequence of these effects is that relatively poor performance can be expected for low flows (ElogQ) when the model is calibrated with an objective function emphasizing high flows (EQ), and vice versa. These results corroborate previous studies dealing with the dependency of the model performance to the type of estimator (i.e. objective function) used during calibration (e.g. Wagener and McIntyre, 2005; Oudin et al., 2006). In particular, the distributed version of mHM parameterized with any of the proposed methods (i.e. HRU and MPR) performed considerably better than the lumped version, regardless of the objective function chosen (Fig. 4, left panels a and c). Furthermore, the performance of the distributed model version exhibited a high sensitivity to the parameterization method employed. The MPR method not only showed less variability but also performed better than the HRU method (Fig. 4, right panels b and d). For instance, the ensemble median of the ElogQ for the distributed model version parameterized with the HRU and the MPR parameterization methods were 12.0% and 18.7% higher than those obtained with the lumped model version during the validation period. Another significant observation drawn from Fig. 4 was that the performance of both models versions obtained with the combined estimator (EQlogQ) was, on average, superior than those obtained with a single objective function (EQ or ElogQ). Moreover, the variability within the good solutions obtained with the combined estimator was, on average, lower than those obtained with the individual estimator for both model versions (Fig. 4). This result was counterintuitive since it was expected that a given model would perform better for high flows if the estimator emphasized only high flows, as compared to those simulations in which the estimator emphasized both high and low flows. A plausible explanation for this effect is related with the water balance and numerical problems such as over-fitting and error propagation. Conditioning a model to match only one kind of event (i.e. high or low flows) would induce bias in many state variables (e.g. soil water content) due to uncertainties in model structure and input. Biased state conditions would, in turn, lead to biased streamflow simulation, and so on.

61

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

(a)

Lumped vs. HRU

HR U v s . M P R

(b) 0.95

0.95

0.85 0.90 0.80

ElogQ

ElogQ

Calibration period

0.90

0.75

0.85 0.70 0.65 0.60

0.80 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95

0.80

0.85

EQ

0.95

0.90

0.95

EQ

(c)

0.95

(d)

0.85

0.90

0.80

ElogQ

Validation period ElogQ

0.90

0.75

0.85

0.70 0.65

0.80

0.60 0.55

0.75 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95

0.75

0.80

EQ

0.90

0.95

EQ

Objective function

0.85

Model Lumped

Dist.- HRU

Dist.- MPR

High flow (EQ) Low flow (ElogQ) Combined (EQlogQ) Fig. 4. Scatterplot of two estimators (EQ and ElogQ) based on daily streamflow simulations at the outlet of study area (i.e. gauge no. 22; Plochingen) for calibration and validation periods. Combinations of lumped and distributed model parameterized with HRU are depicted in panels (a and c), respectively. Comparisons of the distributed model parameterized with HRU and MPR are depicted in panels (b and d), respectively. In all cases, 200 parameters sets per objective function and model version were employed.

Another reason for such performance may be related to the constrained number of model runs per optimization trial, which might have limited the model for achieving a global optimum. It should be noted that for a given model version, the number of model runs per optimization trial was kept constant to have a fair comparison of results. It is worth noting that at least 40,000 model runs were carried out for each model version and objective function. The lumped model version attained the highest level of improvement with the combined estimator (Fig. 4). The distributed version parameterized with the MPR method, on the contrary, improved the least. This result corroborated recent findings (Wagener and McIntyre, 2005; Oudin et al., 2006) regarding the higher sensi-

tivity of lumped models to the type of estimator. Despite this gain the efficiency of the lumped model was, on average, lower than that obtained with the distributed model, regardless of the parameterization method employed. For example, the ensemble median EQ obtained with the distributed model, parameterized with the HRU and the MPR methods, was 2.5% and 7.2% higher than that obtained with the lumped model in the validation period, respectively. This also indicates a clear advantage of the MPR method over the HRU method in this respect. The same tendency in model performance was also observed in other basins as shown in Fig. 5a. The calibration in this case was performed independently for each basin with the combined

62

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

(a)

0.9

EQ

0.8

0.7

median median median P -P

0.6

5

0.5 1

(b)

2

3

95

Lumped Dist.-HRU Dist.-MPR quantile range

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21

0.9 0.8

EQ

0.7 0.6 0.5 0.4 0.3 0.2 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22

1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22

(c) 0.9 0.8

EQ

0.7

0.6

0.5 0.4 0.0 -0.4

Gauge No. Fig. 5. Performance of lumped and distributed (with HRU and MPR parameterizations) mHM for the daily discharge simulations at 22 gauging stations with the split-sample test (panel: a) and the proxy basin test (panels: b and c). Parameter sets obtained through calibration at gauge nos. 22 and 21 were transferred to other locations for the proxy basin test are depicted in panels (b and c), respectively. For all cases calibration were performed with the combined objective function.

estimator as an objective function. The ensemble median of EQ for the distributed model parameterized with MPR was, on average, 15.8% and 5.1% higher than those obtained with the lumped version and the HRU parameterization method, respectively. Similar performance was also observed for the low flow efficiency measure ElogQ. It should be noted that the performance of the lumped model

for various basins considered in this study are consistent with pervious studies reported in the literature (Das et al., 2008; Bárdossy and Singh, 2008). These results also indicated that a critical number of parameters might be required to obtain a good model performance in the investigated basins. The lumped model had the least number of

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

(a)

0.990

r

0.980

0.970

Lump

(b)

Dist.(HRU)

Dist.(MPR)

1.00

to the spatial model discretization and parameterization methods as illustrated in Fig. 6. In this figure, the distributed model, parameterized with the MPR method, performed consistently better than the same model version parameterized with the HRU method. Likewise, the performance of the former was much better than that of the lumped model version. The efficiency measures (r) obtained with the combined estimator, for the lumped and the distributed model versions, were also consistently higher than those obtained with a single estimator. Improvements were particularly remarkable for the cumulative specific deficit as can be seen from Fig. 6b. Hereafter, all sensitivity analysis and remaining results were obtained with the combined estimator. 5.3. Sensitivity of the streamflow simulations to spatial discretization and model parameterization

0.80

r

63

0.60

0.40

0.20 Lump

Dist.(HRU)

Dist.(MPR)

Only high/low flow calibration Combined calibration Ensemble median P5 - P 95 quantile range Fig. 6. Sensitivity of the lumped and the distributed model performances for the prediction of seasonal specific volume of high flows (a) and cumulative specific deficit (b) to the objective functions. The distributed model was parameterized with the HRU and the MPR methods. Simulations were carried out at Plochingen gauging station during the modeling period (1980–2001) with 200 parameter sets found separately for different objective function.

parameters whereas the distributed model version parameterized with the HRU technique had the most. The model with MPR laid in between them. Supporting evidence for this hypothesis was also presented by Das et al. (2008), whom reported that a semi-distributed version of HBV outperformed its lumped and fully distributed versions in the Upper Neckar basin. The number of parameters of the semi-distributed HBV version was also in between that of the distributed and the lumped versions. It is worth noting that the distributed HBV model was parameterized with the HRU method. Further research, however, is still required to test this hypothesis with other models and in other regions. Result obtained with the proxy basin test provided further evidence in favor of the MPR parametrization against that of the HRU method (Fig. 5b and c). The efficiency measures obtained with MPR clearly exhibited less spread (i.e. variance) and higher efficiency (i.e. ensemble median) than those obtained with the HRU method. Furthermore, the MPR method exhibited the least deterioration in performance at presumed ungauged locations as compared to that of the HRU method. This result would imply that the MPR method is more robust and less dependent on calibration location than the HRU method.

5.2. Sensitivity of the model efficiency estimated with extreme runoff characteristics The Pearson correlation coefficient (r) between observed and simulated seasonal specific volume of high flows (Q1) and cumulative specific deficit (Q4) at the outlet of the study area (gauge no. 22) during the modeling period exhibited quite large sensitivity

The ensemble median daily streamflow estimated with the lumped model exhibited large deviations from the observations, specially during recession periods and low flow spells during summer, as illustrated in Fig. 7 for validation water year 1989 at the catchment outlet (gauge no. 22; Plochingen). Moreover, the variability bound (P5–P95) obtained with 200 optimized parameter sets did not cover the streamflow observations in various cases. The parameter calibration for both model versions was performed with the combined objective function. An explanation of this behavior may be related to the lack of spatial variability of basin characteristics in lumped models, which is a major factor controlling recession and the low flow regime in river basins. Due to this reason, the distributed version, regardless of the parameterization employed, exhibited a significant reduction in streamflow variability as compared with the lumped model version. The distributed model parameterized with the MPR method exhibited higher predictive performance as compared to the same model version parameterized with HRU (Fig. 7). A direct result of this was the narrower variability bound (P5–P95) obtained for the simulated streamflow during all recession periods. This effect became even more evident when both parameterization methods were evaluated with the proxy basin test. Fig. 8 shows the daily streamflow simulations for the water year 1988 at gauge no. 21 with model parameters transferred from gauge no. 22. The variability bound of the streamflow simulations with the HRU method was at least 10 % wider than that obtained with the MPR method. Plausible explanations for such large variability in streamflow simulations obtained with the HRU method may be attributed: (1) to the unaccounted sub-grid variability of basin characteristics within an HRU, (2) to the classification method employed to define the HRU classes, (3) to its static character, e.g. land cover changes can not be accounted for, and (4) to the higher number of calibration parameters required by the HRU method as compared to MPR. 5.4. Sensitivity of high flow characteristics to spatial discretization and model parameterization In general, high flow characteristics estimated at the outlet of the study area (gauge no. 22) with the distributed model performed better than the lumped version, regardless of the parameterization method employed (Table 2). The ensemble median RMSE in winter obtained for the distributed model with both HRU and MPR parameterization methods was on average 11% and 24% less than those obtained with the lumped model, respectively. In summer, the same tendency was observed, but the differences were smaller, namely: 8% and 12%, respectively (Table 2). From this table, it can also be noted that the ensemble median RMSE obtained with the HRU method, for

64

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

500

(a)

Observed Ensemble median P5-P95 quantile range

Q(m3s-1)

200

100 50

20

10

500

(b)

Q(m3s-1)

200

100 50

20

10

500

(c)

Q(m3s-1)

200

100 50

20

10

Nov

Dec

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Water year 1989 Fig. 7. Ensemble of daily streamflow simulations at the outlet of study area (i.e. gauge no. 22; Plochingen) during the validation water year 1989 using the lumped model (a) and the distributed model with the HRU (b) and HRU (b) and the MPR (c) parameterization methods. The variability range (P5–P95) of the modeled streamflow were estimated with 200 parameter sets. These sets were obtained with the combined objective function.

300 200

(a)

Observed Ensemble median P5-P95 quantile range

Q(m3s-1)

100 50 30 20

10 5 3 2 300 200

(b)

Q(m3s-1)

100 50 30 20

10 5 3 2

Nov

Dec

Jan

Feb

Mar

Apr

May

Jun

Jul

Aug

Sep

Oct

Nov

Water year 1988 Fig. 8. Ensemble of daily streamflow simulations at Horb gauging station (gauge no. 21) during the water year 1988 using the HRU (a) and the MPR (b) parameterization methods. The ensemble median and variability range (P5–P95) of the modeled streamflow were estimated with the 200 parameter sets. These sets were obtained with the combined objective function and the streamflow data of Plochingen gauging station (gauge no. 22).

65

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

Table 2 Ensemble median statistics between observed and simulated seasonal high flow characteristics at Plochingen (gauge no. 22) and Horb (gauge no. 21) during the modeling period (1980–2001). The lumped model and the distributed model (HRU and MPR parameterizations) simulations were obtained with 200 parameter sets. Simulations corresponding to Horb were based on the proxy basin test with parameters sets transferred from calibration at Plochingen. Season/characteristic

Plochingen

Horb

RMSE

r

RMSE

r

Lumped

Dist. HRU

Dist. MPR

Lumped

Dist. HRU

Dist. MPR

Dist. HRU

Dist. MPR

Dist. HRU

Dist. MPR

Winter Q1 (mm) Q2 (d) Q3 (y1)

13.34 3.49 1.37

12.11 3.04 1.21

9.93 2.61 1.07

0.97 0.95 0.85

0.98 0.96 0.89

0.99 0.98 0.91

29.73 6.69 1.65

20.94 4.65 1.07

0.95 0.94 0.83

0.97 0.96 0.87

Summer Q1 (mm) Q2 (d) Q3 (y1)

9.01 2.01 1.05

8.05 1.93 0.96

7.96 1.87 0.86

0.90 0.89 0.86

0.90 0.89 0.88

0.91 0.90 0.89

7.80 1.38 0.83

7.03 1.11 0.69

0.85 0.84 0.82

0.88 0.90 0.90

winter and summer, was on average 17% and 5% greater than those obtained with the MPR method, respectively. However, the Pearson’s correlation (r) obtained with the distributed model, for winter and summer respectively, did not exhibit large differences as compared with those obtained with the lumped model. The largest difference found (i.e. 7%) corresponded to the r value between the distributed model using MPR and the lumped version for the frequency of high flows (Q3). Lumped and distributed model versions showed relatively poor performance for the predictions of the frequency of high flows (Q3) as compared to the specific volume (Q1) and the total duration (Q2) of high flows for both seasons. In winter, the r value was on average 10% less than those obtained for the other two runoff characteristic presented in Table 2. In summer, on the contrary, this difference was 2% smaller. It should be noted, however, that this relative reduction was mainly due to the poor performance of the specific volume and the total duration in summer. Additionally, this behavior was observed regardless of the parameterization method employed for the distributed model. This effect can also be visually observed from Fig. 9. A significant difference in performance was observed between lumped and distributed model versions when they were evaluated

at remaining locations. Fig. 10 depicts the mean, maximum and minimum range of the ensemble median correlation coefficient r obtained for the remaining basins. The correlation coefficient r between the observed and the simulated seasonal high flow characteristics obtained with the distributed model clearly exhibited less spread and higher efficiency than those obtained with the lumped model. Moreover, the distributed model parameterization based on the MPR method showed, on average, better performance than that obtained with the HRU method. Model performances were generally lower during summer than those in winter. This behavior was noticed regardless of the type of spatial discretization or parameterization method employed (Fig. 10). It is worth noting that similar results have been also obtained with parsimonious nonlinear generalized models for the same study area (Samaniego and Bárdossy, 2007). The lack of predictability of high flow characteristics in summer is likely related with the occurrence of convective precipitation events, whose significant large intensity, relatively short duration, and limited spatial extent is not often captured by the network of weather stations available in this region. The poor performance of the lumped model version compared to those of the distributed version during summer also supported this hypothesis. The performance of the

Observed

12

12

10

10

10

8

)

8

)

8

6

6

Q (y

6

4

4

0

0

50

50

50

40

40

40

20

30 20

30 20

10

0

0

0

250

250

250

200

200

200

100

10

Q (mm)

10

150

P5 - P95 quantile range

2

Q (d)

0

Q (d)

2

30

Ensemble median

4

2

Q (mm)

Q (mm)

Q (d)

12

Q (y

(c)

)

(b)

Q (y

(a)

150 100

150 100

50

50

50

0

0

0

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

Water year

Water year

Water year

Fig. 9. Performance of the lumped (a) and the distributed model with the HRU (b) and the MPR (c) parameterization methods for the prediction of seasonal high flow characteristics at the outlet of study area (i.e. gauge no. 22; Plochingen) in the modeling period (1980–2001).

66

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

The robustness of the MPR method over the HRU method was also supported by the proxy basin test carried out at Horb (gauge no. 21) using the parameters sets obtained at Plochingen (gauge no. 22) (Table 2). For instance, the median RMSE values between observed and simulated high flow characteristics obtained with the MPR were on average 45% and 18% less than those values obtained with the HRU method for winter and summer, respectively. The ensemble median r also showed the similar tendency.

1

0.8

r

0.6

0.4

mean Lumped mean Dist. - HRU mean Dist. - MPR

0.2

5.5. Sensitivity of low flow characteristics to spatial discretization and model parameterization

max - min range

0 Q1

Q2

Q3

Q1

Q2

Winter

Q3

Summer

Seasonal high flow characteristics Fig. 10. Mean and range of ensemble median r obtained between observed and simulated seasonal high flow characteristics using lumped and the distributed (HRU and MPR) model versions. Both statistics were obtained from the model runs of 21 basins (except gauge no. 22) for the modeling period (1980–2001).

extreme high flow characteristics obtained with mHM were, nevertheless, better than those obtained with simple data-based driven models as reported by Samaniego and Bárdossy (2007).

Low flow characteristics estimated at the outlet of the study area (gauge no. 22) with the distributed model also performed better than the lumped version, regardless of the parameterization method employed (Table 3). The ensemble median RMSE in winter, obtained for the distributed model with both HRU and MPR parameterization methods, was on average 29% and 38% less than those obtained with the lumped model, respectively. In summer, the same tendency was observed, namely: 20% and 41%, respectively (Table 3). From this table, it can also be noted that the ensemble median RMSE obtained with the HRU method, for winter and summer, was on average 12% and 18% greater than those obtained with the MPR method, respectively.

Table 3 Ensemble median statistics between observed and simulated seasonal low flow characteristics at Plochingen (gauge no. 22) and Horb (gauge no. 21) during the modeling period (1980–2001). The lumped model and the distributed model (HRU and MPR parameterizations) simulations were obtained with 200 parameter sets. Simulations corresponding to Horb were based on the proxy basin test with parameters sets transferred from calibration at Plochingen. Season/characteristic

Plochingen

Horb

RMSE

r

RMSE

Lumped

Dist. HRU

Dist. MPR

Lumped

Dist. HRU

Dist. MPR

Winter Q4 (mm) Q5 (d) Q6 (mm y1)

0.58 8.14 9.69

0.41 5.57 7.16

0.32 5.03 6.51

0.80 0.76 0.78

0.82 0.85 0.87

0.86 0.90 0.93

Summer Q4 (mm) Q5 (d) Q6 (mm y1)

0.56 8.51 7.03

0.45 7.62 4.21

0.37 4.64 3.87

0.82 0.81 0.69

0.93 0.92 0.88

0.96 0.94 0.92

(a)

(b)

60

Dist. HRU

Dist. MPR

Dist. HRU

Dist. MPR

0.41 9.06 6.23

0.23 5.21 5.70

0.84 0.88 0.92

0.96 0.95 0.97

0.92 11.49 8.27

0.63 6.58 4.13

0.89 0.86 0.83

0.93 0.92 0.95

(c)

60

r

60

Observed Ensemble median

20

40

Q (mm y

Q ( mm y

40

20

20

0 80

60

60

60

40

Q (d)

0 80

Q (d)

0 80

40

40

20

20

0

0

0

5

5

5

4

4

4

3 2

Q ( mm)

20

Q (mm)

Q 6( mm y Q (d) Q (mm)

)

)

)

P5 - P95 quantile range 40

3 2

3 2

1

1

1

0

0

0

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

Water year

Water year

Water year

Fig. 11. Performance of the lumped (a) and the distributed model with HRU (b) and MPR (c) parameterization methods for the prediction of seasonal low flow characteristics at the outlet of study area (i.e. gauge no. 22; Plochingen) in the modeling period (1980–2001).

67

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

with the lumped model, namely 10% and 19%, respectively. The MPR method in this case also outperformed the HRU method in winter and summer. Predictions made with the lumped model version not only significantly mismatched observations but also exhibited the largest variability (P5–P95) as compared to those obtained with the distributed model version (Fig. 11). Similar behavior was also noticed at other locations as depicted in Fig. 12. The distributed model version, regardless of the parameterization method, performed, on average, better than the lumped version, specially during summer. Differences in model performance due to the parameterization methods were also apparent. The distributed model parameterization with the MPR method in this case also exhibited better performance than that obtained with the HRU method. The results of the proxy basin test carried out at Horb (gauge no. 21) with parameter sets calibrated with the data of Plochingen (gauge no. 22) further corroborated the robustness and reliability of MPR over HRU (Table 3). For instance, the ensemble median RMSE obtained with the MPR was, on average, 32% lower than that obtained with the HRU method. The variability bound of the ensemble low flow characteristics obtained with MPR was also significantly narrower than that obtained with the HRU method as depicted in Fig. 13. The relatively poor performance of the lumped model as compared to spatially distributed model versions (Figs. 11 and 12) highlighted the effects of the spatial discretization and thus the importance of accounting the spatial heterogeneity of catchment characteristics in modeling of low flow characteristics in the investigated basins. Results supporting this hypothesis could also be inferred by the better performance of the MPR method as compared to the HRU method. The MPR based parametrization method accounts for the sub-grid variability of the catchment characteristics in the first step of the regionalization procedure, whereas the HRU method uses the grid scale (i.e. at modeling scale, level-2) catchment characteristics in the HRUs classification scheme. It may be

Moreover, the RMSE value of the total drought duration in summer (Q5) obtained with the lumped model version is only 9% less than that obtained with the best input–output relationship found by Samaniego and Bárdossy (2007) for the same study area. The distributed model using MPR, on the other hand, produced a reduction of almost 100% as compared to the same benchmark. The Pearson’s correlation (r) obtained with the distributed model, regardless of the parameterization method employed, was for both winter and summer, significantly better than that obtained

mean Lumped mean Dist. - HRU mean Dist. - MPR max - min range 1 0.8

r

0.6 0.4 0.2 0 Q4

Q5

Q6

Q4

Winter

Q5

Q6

Summer

Seasonal low flow characteristics Fig. 12. Mean and range of ensemble median r obtained between observed and simulated seasonal low flow characteristics using lumped and the distributed (HRU and MPR) model versions. Both statistics were obtained from the model runs of 21 basins (except gauge no. 22) for the modeling period (1980–2001).

60

60

40 20

20 0

100

80

80

Q5(d)

0

60 40

60 40

20

20

0 10

0 10

8

8

6 4 2

Observed Ensemble median P5-P95 quantile range

40

100

Q4(mm)

Q5(d) Q4(mm)

Q6(mm y -1)

(b) 80

Q6(mm y -1)

(a) 80

6 4 2

0

0 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002

Water year

Water year

Fig. 13. Performance of the HRU (a) and the MPR (b) parameterization methods employed in the distributed mHM model for the prediction of seasonal low flow characteristics at Horb in the modeling period (1980–2001).

68

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69

noted that the differences in distributed model performance can only be attributed to the changes in the parameterization method since the model structure and the meteorological forcing were kept same for both parameterization methods.

6. Conclusions The purpose of this study was to investigate the effects of the spatial discretization, parameterization method, and type of objective function, not only on the daily streamflow simulations but also on seasonal extreme runoff characteristics. To address this goal, lumped and distributed versions of mHM were applied. The distributed version was parameterized, in turn, with two methods: HRU and MPR. Good sets of model parameter were identified through calibration with three objective functions emphasizing: high flows, low flows and the combination of both. The following conclusions were drawn based on results obtained for 22 German basins. The performance of mHM for the daily streamflow simulations was sensitive to the type of estimator employed during calibration, regardless of the spatial discretization and the parameterization method. Both model versions calibrated with the combined estimator exhibited higher model efficiency and less spread than those obtained with a single estimator. The model efficiency was also sensitive to the spatial discretization and parameterization method. The distributed model version parameterized with the MPR method consistently attained the best performance irrespective of the estimator employed, whereas the lumped model version always attained the least performance. The model efficiency based on seasonal runoff characteristics also exhibited a strong sensitivity to spatial discretization and model parameterization. Likewise, the distributed model parameterized with MPR and the lumped model version attained the best and the worst performance, respectively. Numerical modeling results provided compelling evidence supporting the role of the spatial variability and sub-grid heterogeneity of catchment characteristics on recession periods and low flow spells during summer. In these cases, the distributed model version parameterized with the MPR method consistently outperformed the distributed model parameterized with the HRU method as well as the lumped version. All six extreme runoff characteristics investigated in this study exhibited marked sensitivity to the spatial discretization and the parameterization method. In this respect, the distributed model version consistently outperformed the lumped model, regardless of the parameterization method used. A significant improvement in the predictability of both high and low flow characteristics was noticed when the distributed model was used instead of the lumped version, which clearly highlighted the advantages of distributed modeling. The MPR method not only reduced the number of free calibration parameters but also exhibited significant improvement in the predictability of extreme runoff characteristics as compared to those obtained through the HRU method. The proxy basin test consistently supported the MPR method with respect to robustness and reliability for the simulation of daily streamflow and extreme runoff characteristics over the HRU method. This result may be relevant for predictions in ungauged basins (PUB). Acknowledgments The authors would like to thank the Editor, the Associate Editor, Dr. Jens Götzinger and four anonymous reviewers for their valuable suggestions and critique on the manuscript.

References Ahmed, S., Marsily, G.d., 1987. Comparison of geostatistical methods for estimating transmissivity using data on transmissivity and specific capacity. Water Resour. Res. 23 (9), 1717–1737. Bárdossy, A., Singh, S.K., 2008. Robust estimation of hydrological model parameters. Hydrol. Earth Syst. Sci. 12 (6), 1273–1283. Beldring, S., Engeland, K., Roald, L.A., Saelthun, N.R., Vokso, A., 2003. Estimation of parameters in a distributed precipitation–runoff model for Norway. Hydrol. Earth Syst. Sci. 7 (3), 304–316. Bergström, S., 1995. The HBV model. In: Singh, V. (Ed.), Computer Models of Watershed Hydrology. Water Resour. Publ., Colorado, USA, pp. 443–476. Beven, K., 2006. A manifesto for the equifinality thesis. J. Hydrol. 320 (1–2), 18–36. Blöschl, G., Reszler, C., Komma, J., 2008. A spatially distributed flash flood forecasting model. Environ. Model. Softw. 23 (4), 464–478. Booij, M.J., 2005. Impact of climate change on river flooding assessed with different spatial model resolutions. J. Hydrol. 303, 176–198. Bronstert, A., Niehoff, D., Bürger, G., 2002. Effects of climate and land-use change on storm runoff generation: present knowledge and modelling capabilities. Hydrol. Process. 16 (2), 509–529. Das, T., Bárdossy, A., Zehe, E., He, Y., 2008. Comparison of conceptual model performance using different representations of spatial variability. J. Hydrol. 356, 106–118. Davison, B., van der Kamp, G., 2008. Low flows in deterministic modeling: a brief review. Can. Water Resour. J. 33 (2), 181–194. Dickinson, R., 1984. Modelling evapotranspiration for three-dimensional global climate models. In: Hansen, J.E., Takahashi, T. (Eds.), Climate Processes and Climate Sensitivity, Geophysical Monograph Series, vol. 29. AGU, pp. 58–72. Flügel, W.A., 1995. Delineating hydrological response units by geographical information system analyses for regional hydrological modelling using PRMS/ MMS in the drainage basin of the river Bröl, Germany. Hydrol. Process. 9, 423– 436. Götzinger, J., Bárdossy, A., 2007. Comparison of four regionalisation methods for a distributed hydrological model. J. Hydrol. 333, 374–384. Grayson, R., Blöschl, G., 2000. Spatial Patterns in Catchment Hydrology: Observations and Modelling. Cambridge University Press, The Edinburgh Building, Cambridge, UK. ISBN 0-521-63316-8. Gupta, H.V., Sorooshian, S., Yapo, P., 1998. Toward improved calibration of hydrologic models: multiple and noncommensurable measures of information. Water Resour. Res. 34 (4), 751–763. Hargreaves, G.H., Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Appl. Eng. Agri. 1, 96–99. Hundecha, Y., Bárdossy, A., 2004. Modeling effect of land use changes on runoff generation of a river basin through parameter regionalization of a watershed model. J. Hydrol. 292, 281–295. Koren, V., Finnerty, R., Schaake, J., Smith, M., Seo, D., Duan, Q., 1999a. Scale dependencies of hydrologic models to spatial variability of precipitation. J. Hydrol. 217, 285–302. Koren, V., Schaake, J., Mitchell, K., Duan, Q.Y., Chen, F., Baker, J.M., 1999b. A parameterization of snowpack and frozen ground intended for NCEP weather and climate models. J. Geophys. Res. 104, 19569–19585. Krajewski, W., Lakshmi, V., Georgakakos, K., Jain, S., 1991. A Monte Carlo study of rainfall sampling effect on a distributed catchment model. Water Resour. Res. 27 (1), 119128. Kundzewicz, Z.W., Mata, L.J., Arnell, N.W., Döll, P., Kabat, P., Jiménez, B., Miller, K.A., Oki, T., Sen, Z., Shiklomanov, I.A., 2007. Freshwater resources and their management. In: Parry, M.L., Conziani, O.F., Palutikof, J.P., Van der Linden, P.J., Hanson, C.E. (Eds.), Climate Change 2007: Impacts, Adaptation and Vulnerability, Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on climate Change. Cambridge University Press, Cambridge, UK. Kutilek, M., Nielsen, D., 1994. Soil Hydrology. Catena Verlag, Cremlingen-Destedt. LeMoine, N., Andréassian, V., Perrin, C., Michel, C., 2007. How can rainfall–runoff models handle intercatchment groundwater flows? Theoretical study based on 1040 french catchments. Water Resour. Res. 43, W06428. Leavesley, G.H., Lichty, R.W., Troutman, B.M., Saindon, L.G., 1983. Precipitation– runoff Modeling System: User’s Manual, US Geological Survey Water-Resources Investigations. Denver – Colorado, Rep. 83–4238. Lehner, B., Döll, P., Alcamo, J., Henrichs, H., Kaspar, F., 2006. Estimating the impact of global change on flood and drought risks in Europe: a continental integrated analysis. Clim. Change 75, 273–299. Liang, X., Lettenmaier, D.P., Wood, E.F., Burges, S.J., 1994. A simple hydrologically based model of land-surface water and energy fluxes for general-circulation models. J. Geophys. Res. 99 (D7), 14415–14428. Middelkoop, H., Daamen, K., Gellens, D., Grabs, W., Kwadijk, J., Lang, H., Parmet, B., Schädler, B., Schulla, J., Wilke, K., 2001. Impact of climate change on hydrological regimes and water resources management in the Rhine Basin. Clim. Change 49 (1), 105–128. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models, 1. A discussion of principles. J. Hydrol. 10, 282–290. Ouarda, T., Charron, C., St. Hilaire, A., 2008. Statistical models and the estimation of low flows. Can. Water Resour. J. 33 (2), 195–206. Oudin, L., Andréassian, V., Perrin, C., Michel, C., 2006. Dynamic averaging of rainfall– runoff model simulations from complementary model parameterizations. Water Resour. Res. 42, W07410.

R. Kumar et al. / Journal of Hydrology 392 (2010) 54–69 Patterson, D.E., Smith, M.W., 1981. The measurement of unfrozen water content by time domain reflectometry: results from laboratory tests. Can. Geotech. J. 18, 131–144. 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, 275–301. Pilgrim, D.H. (Ed.), 1987. Australian Rainfall and Runoff, a Guide to Flood Estimation. Inst. of Eng. Aust., Canberra. Pokhrel, P., Gupta, H.V., Wagener, T., 2008. A spatial regularization approach to parameter estimation for a distributed watershed model. Water Resour. Res. 44, W12419. Reed, S., Koren, V., Smith, M., Zhang, Z., Moreda, F., Seo, D., Participants, D., 2004. Overall distributed model intercomparison project results. J. Hydrol. 298 (1–4), 27–60. Samaniego, L., 2003. Hydrological Consequences of Land Use/Land Cover Change in Mesoscale Catchments. Institute of Hydraulic Engineering, University of Stuttgart, Faculty of Civil Engineering, Stuttgart, Ph.D. Dissertation No. 118, ISBN 3-9337 61-21-2. Samaniego, L., Bárdossy, A., 2007. Relating macroclimatic circulation patterns with characteristics of floods and droughts at the mesoscale. J. Hydrol. 335, 109–123. Samaniego, L., Bárdossy, A., Kumar, R., 2010a. Streamflow prediction in ungauged catchments using copula-based dissimilarity measures. Water Resour. Res., W02506. Samaniego, L., Kumar, R., Attinger, S., 2010b. Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale. Water Resour. Res., W05523.

69

Seibert, J., 2003. Reliability of model predictions outside calibration conditions. Nord. Hydrol. 34, 477–492. Shevenell, L., 1999. Regional potential evapotranspiration in arid climates based on temperature, topography and calculated solar radiation. Hydrol. Process. 13 (13), 577–596. Smakhtin, V.U., 2001. Low flow hydrology: a review. J. Hydrol. 240, 147–186. Stedinger, J.R., Vogel, R.M., Foufoula-Georgiou, E., 1993. Frequency analysis of extreme events. In: Handbook of Hydrology. McGraw-Hill, New York, USA. Tewolde, M.H., Smithers, J.C., 2006. Flood routing in ungauged catchments using Muskingum methods. Water SA 32 (3), 379–388. Tolson, B.A., Shoemaker, C.A., 2007. Dynamically dimensioned search algorithm for computationally efficient watershed model calibration. Water Resour. Res. 43, WR004723. Tolson, B.A., Shoemaker, C.A., 2008. Efficient prediction uncertainty approximation in the calibration of environmental simulation models. Water Resour. Res. 44, WR005869. Uhlenbrook, S., Roser, S., Tilch, N., 2004. Hydrological process representation at the meso-scale: the potential of a distributed, conceptual catchment model. J. Hydrol. 291 (3–4), 278–296. Wagener, T., McIntyre, N., 2005. Identification of rainfall–runoff models for operational applications. Hydrol. Sci. J. 50 (5), 1–18. Zacharias, S., Wessolek, G., 2007. Excluding organic matter content from pedotransfer predictors of soil water retention. Soil. Sci. Soc. Am. J. 71 (1), 43–50.