Watershed model parameter estimation and uncertainty in data-limited environments

Environmental Modelling & Software 51 (2014) 84e93

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

Watershed model parameter estimation and uncertainty in data-limited environments André Fonseca a, *, Daniel P. Ames b, Ping Yang c, Cidália Botelho a, Rui Boaventura a, Vítor Vilar a, * a LSRE e Laboratory of Separation and Reaction Engineering e Associate Laboratory e LSRE/LCM, Department of Chemical Engineering, Faculty of Engineering of the University of Porto, Porto, Portugal b Department of Civil and Environmental Engineering, Brigham Young University, Provo, UT, USA c Texas Institute for Applied Environmental Research, Tarleton State University, Stephenville, TX, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 3 June 2013 Received in revised form 25 September 2013 Accepted 27 September 2013 Available online

Parameter uncertainty and sensitivity for a watershed-scale simulation model in Portugal were explored to identify the most critical model parameters in terms of model calibration and prediction. The research is intended to help provide guidance regarding allocation of limited data collection and model parameterization resources for modelers working in any data and resource limited environment. The watershed-scale hydrology and water quality simulation model, Hydrologic Simulation Program e FORTRAN (HSPF), was used to predict the hydrology of Lis River basin in Portugal. The model was calibrated for a 5-year period 1985e1989 and validated for a 4-year period 2003e2006. Agreement between simulated and observed streamflow data was satisfactory considering the performance measures such as NasheSutcliffe efficiency (E), deviation runoff (Dv) and coefficient of determination (R2). The Generalized Likelihood Uncertainty Estimation (GLUE) method was used to establish uncertainty bounds for the simulated flow using the NasheSutcliffe coefficient as a performance likelihood measure. Sensitivity analysis results indicate that runoff estimations are most sensitive to parameters related to climate conditions, soil and land use. These results state that even though climate conditions are generally most significant in water balance modeling, attention should also focus on land use characteristics as well. Specifically with respect to HSPF, the two most sensitive parameters, INFILT and LZSN, are both directly dependent on soil and land use characteristics. Ó 2013 Elsevier Ltd. All rights reserved.

Keywords: HSPF GLUE Monte Carlo simulations Watershed modeling Uncertainty assessment Sensitivity analysis

Software and data availability Name

BASINS 4.0 (Better Assessment Science Integrating point & Non-point Sources) with a non-proprietary, open source, free GIS system, MapWindow (www.MapWindow.org) Developer U.S. EPA with AquaTerra Consultants and Idaho State University Contact http://www.aquaterra.com/contact/index.php Availability and cost The software is available for free download at USEPA (United States Environmental Protection Agency) website. Mapwindow is an open source programmable GIS (VB, Cþþ, .NET, and Active X controls) that supports manipulation, analysis, and viewing of geospatial data and associated attribute data in several GIS data formats. Name GLUEWIN Developer EEMC e Euro-area Economy Modelling Centre Contact [email protected] Availability and Cost GLUEWIN is a code designed for analyzing the output of Monte Carlo runs when empirical observations of the model output are available and implements the combination of GSA and GLUE methodologies. The software can be obtained for free at EEMC website (http://eemc.jrc.ec.europa.eu/Software-GLUEWIN.htm).

* Corresponding authors. Tel.: þ351 918257824; fax: þ351 225081674. E-mail addresses: [email protected] (A. Fonseca), [email protected] (V. Vilar). 1364-8152/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.envsoft.2013.09.023

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

1. Introduction Like many similar government initiatives throughout the world, the European Union Water Framework Directive (WFD) was established to restore and protect both surface and ground water with ambitious goals to be met by a target date of 2015 (2000/60/ EC, 2000). Watershed modeling software can be used to help scientists and watershed managers to meet these goals by simulating the effect on water quality and quantity considering different water management strategies, types of land use and climate change. Many commercial and open source watershed simulation models are available and much care and consideration needs to be employed when choosing a model for application in a particular watershed (Borah and Bera, 2004). Regardless of the model chosen, even greater attention must be given to the task of “populating” the selected model with appropriate and physically meaningful model parameters that accurately characterize the surficial landscape, subsurface geology, atmospheric conditions, and other constraints affecting the storage and flux of water and contaminants through the environment (Donigian, 2002; Doherty and Johnston, 2003). Because many watershed simulation model parameters are difficult or impossible to measure in the natural world, parameters must often be estimated or otherwise evaluated from secondary information sources and hence are typically laden with notable degrees of uncertainty (Gallagher and Doherty, 2006). This paper presents a watershed modeling study in the Lis River basin (Portugal) with the express purpose of identifying those parameters in the selected watershed model, whose accurate characterization is most critical for successful model application. The study includes a complete model parameterization and calibration effort combined with parameter uncertainty estimation techniques. Great attention is given to accurate characterization of those parameters for the success of the modeling effort. The results provided here can be used to inform other modelers in data and resource limited situations as to which parameters warrant the greatest resource-allocation and technical attention e allowing for the more efficient use of limited resources. The hydrological model, Hydrologic Simulation Program FORTRAN (HSPF) was used in this study. HSPF is based on the original Stanford Watershed Model IV (Crawford and Linsley, 1966) and is a consolidation of three previously developed models: Agricultural Runoff Management Model (ARM) (Donigian and Davis, 1978), Non-point Source Runoff Model (NPS) (Donigian and Crawford, 1976) and Hydrological Simulation Program (HSP) including HSP Quality (Donigian et al., 1991, 1995; Hydrocomp, 1977). HSPF is a semi-distributed model that simulates water and contaminant transport through spatially distributed, physically homogenous areas within a watershed called Hydrologic Response Units (HRUs). HRUs are presumed to hydrologically respond similarly to given meteorological inputs (precipitation, potential evapotranspiration and temperature). In this way, HSPF can simulate the hydrological, hydraulic and water quality processes on pervious and impervious land surfaces, in soil profiles and in streams and well-mixed impoundments on a continuous basis (Bicknell et al., 2001). A graphical user interface for HSPF is included in the free software, BASINS, developed and distributed by the United States Environmental Protection Agency (http://water.epa.gov/scitech/ datait/models/basins/index.cfm). BASINS is built on the open source geographic information system (GIS) MapWindow GIS (Ames et al., 2008). The use of BASINS to develop HSPF models has been reported in several studies (Bergman et al., 2002; Carrubba,

85

2000; Lian et al., 2007; Lowe and Doscher, 2003; J. Zhang et al., 2009). Calibration of HSPF is an iterative procedure of parameter evaluation, as a result of comparing simulated against observed values of interest. Since HSPF uses a large number of parameters that can be adjusted to represent the physical environment, an expert system, HSPExp is available to assist modelers with hydrology calibration. Typically a dozen or less parameters are used in most studies. HSPExp advises the user on which parameters can be meaningfully adjusted to reduce simulation error while providing explanations regarding the modifications (Lumb et al., 1994). However, there are limitations to HSPF, such as limited spatial definition (finite element analysis model), limited to non-tidal freshwater systems and extensive data requirements (i.e. meteorological and most important gaging stations of interest in the watershed). Model validation is necessary in any model application and is an extension of the calibration process. Its purpose is to assure that the calibrated model properly assesses all the variables and conditions which can affect model results and, the ability to predict the behavior for periods separate from the calibration. Model credibility is based on the ability of a single set of parameters to represent the entire range of observed data. If a single parameter set can reasonably represent wide range of events, then this is a form of validation (Donigian, 2002). Sensitivity and uncertainty analysis of model parameters is conventionally considered to be one of the primary steps in the development and evaluation of models (Sudheer et al., 2011; Jakeman et al., 2006). Over the past decade it has become widely accepted that hydrological models with numerous parameters are likely to produce equally acceptable predictions for multiple different parameter sets (Hope et al., 2004) and a unique “best” parameter set cannot necessarily be found in the parameter space (Christiaens and Feyen, 2001). A structured method to quantify model uncertainty, Generalized Likelihood Uncertainty Estimation (GLUE), proposed by K. Beven and Binley (1992) establishes uncertainty bounds for the simulated value using parameter sets that are determined to be acceptable based on a performance likelihood measure. The implementation of this method requires the user to make a number of subjective decisions such as the threshold value of the likelihood measure to classify the model as acceptable or unacceptable (K. Beven and Binley, 1992; K. J. Beven, 2001). With GLUE, model parameters are sampled from distributions, typically with independent uniform or normal distributions for each parameter. The model is then run with each parameter set, generating multiple sets of model output, which are used to generate uncertainty intervals for model predictions. The generated model parameters are grouped in two categories: behavioral, sets of model parameters that produce results consistent with the observations and non-behavioral, results that contradict the observations (Spear and Hornberger, 1980). As discussed previously the model HSPF, requires extensive parameterization to simulate hydrology fluxes due to both uncertainty in modeled processes and observation errors (Ratto et al., 2001). The remainder of this paper includes a discussion of the methods employed, and results related to the hydrology calibration parameters with focus on assessing the sensitivity of the model to selected parameters and determination of parameter priority in the model. 2. Methods 2.1. Study area The Lis River basin is one of the most important natural resources of the Leiria region in Portugal, with a population rich in fish species and much sought for sport

86

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

fishing events. It covers an area of approximately 853 km2 and is located in the central region of Portugal. It is a coastal basin that is confined to the north by the basin of the Mondego River, and to the south by the Alcoa River basin (Fig. 1). Meteorological and discharge data for the Lis River basin was retrieved from SNIRH e Sistema Nacional de Informação de Recursos Hídricos (National Information System for Water Resources). A high resolution Digital Elevation Model (DEM, 25 m) was provided by ARH-Centro, Administração de Recursos Hidricos (Water Resources Administration e Center Portugal) and the Land Use and Soil Data from Agência Portuguesa do Ambiente- Atlas do Ambiente (Ministry of Agriculture, Sea, Environment and Spatial Planning). Table 1 shows the climate and geographic characteristics of the basin. According to hydrologic soil groups classification (Donigian and Davis, 1978) it is considered to have a moderate and moderate to high runoff potential (Group B and C). Land use categories were aggregated into 5 main categories (Mapoteca, 2008), to simplify the modeling process (Urban or Built-up Area, Forest Land, Agricultural Land, Barren Land and Wetlands) (Table 2). 2.2. HSPExp modeling approach HSPF, coupled with BASINS data management and graphical user interface tools, was used to model the hydrology of the Lis River watershed in Leiria (Portugal) for the years 1985e1989. As noted previously, HSPExp provides the user with expert advice on which parameters should be changed up or down to improve the calibration (Lumb et al., 1994). A list of hydrology and hydraulic parameters and corresponding recommended value ranges can be found in BASINS Technical Note 6 (EPA, 2000). HSPExp uses over 35 rules involving over 80 conditions to recommend parameter adjustments. These rules are divided into four phases: Water balance; High/Low flow; Storm flow; Seasonal adjustments (Lumb et al., 1994). After one phase meets the criteria established, HSPExp moves to the next phase and so on until all criteria are met. 2.3. Simulation data and efficiency assessment The precipitation dataset used was interpolated from 2 stations located within the watershed using the BASINS WDMUtil software tool and were converted from daily to hourly time step from 1985 to 1989 as required by the HSPF platform. The simulation was performed at the gaging station 15E/03, Lena River with a catchment area of 176 km2 (Fig. 2), a main tributary of Lis River where observed data is available for the 5 year period. Lena River basin is comprised of 46% Agriculture, 24% Forest, 20% Barren and 10% Urban (wetland omitted, approximately 0.6%). The calibration approach was independent of land use characteristics except for the index to lower zone evapotranspiration (LZETP) which according to the BASINS Technical Note 6 different values are expected for each land use. The flow data available at station 15E/03 is measured as hydrometric level, and only one flow curve was available down to a minimum hydrometric level of 0.60379. Some of the observed data was below that level where we could not calculate flow. For those days, period of year 1988, we simply considered zero flow rate. A total of 18 storms were selected by observing the plots of daily precipitation and daily mean flow. HSPExp allows entry up to 36 major storms, which selection should be distributed evenly throughout the year and the storm period should start on the first day of significant precipitation and extend for several days to the time when flow returns to near pre-storm levels and/or the hydrograph flattens out

Table 1 Climate conditions and geographic characteristics of Lis Basin. Climate conditions Annual average precipitation (mm) Annual average temperature ( C) Snow is not a factor

Geographic characteristics 855

Average slope (%)

14.8

9.6

Soil (%): Silt loam Sandy clay loam Sandy loam

46 40 14

significantly (USEPA). Daily potential evapotranspiration was derived with the WDMUtil tool using maximum and minimum daily air temperature values since it was not available at any meteorological station in Lis watershed. The average percentage for impervious amount of area was defined according to the values given in the land use shapefile (Table 3). After an HSPF project is created, hydrology (water quantity) is generally calibrated before any water quality calibration effort. The performance of the model for simulating runoff was evaluated using the HSPExp criteria and statistical methods, the NasheSutcliffe coefficient of efficiency (E), the Deviation Runoff (Dv) and the Coefficient of Determination (R2). The NasheSutcliffe coefficient of efficiency (Nash and Sutcliffe, 1970), E, is a widely used and potentially reliable statistic for assessing the goodness of fit of hydrologic models, and it is defined as: E ¼ 1


Xn i¼1

ðOi  Si Þ2

.
Xn i¼1

ðOi  BÞ2



where O is the measure streamflow, S is the simulated streamflow, Ø is the average measure streamflow and n the number of streamflow values. E can range from N to 1. An efficiency of 1 (E ¼ 1) corresponds to a perfect match between observed and predicted data, and values equal or less than zero indicate that the model is predicting no better than using the average of the observed data (ASCE, 1993). E is frequently used in the hydrological literature for the evaluation of model performance (Baloch and Tanik, 2001; Jin et al., 2010; Johnson et al., 2003; Li et al., 2010; Zhang et al., 2009). The Deviation of Runoff (Dv) calculation allows the absolute differences to be weighted by the actual runoff volumes, with a value of zero indicating an ideal model fit. The Dv is given by the following equation (Martinec and Rango, 1989): Dv ½% ¼

h
Xn i¼1

.
Xn i ðOi  Si Þ ðOi Þ  100 i¼1

where O is the measure streamflow, S is the simulated streamflow and n the number of streamflow values. Dunker and Melching (1998) suggest that HSPF model calibration can be considered satisfactory if the E value for monthly flows exceeds 0.80. According to Donigian et al. (1984), monthly and annual simulations are very good when percentage error (Dv) is less than 10, good when it is between 10 and 15, and fair when it is between 15 and 25.

Fig. 1. Lis watershed in Leiria, Portugal.

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93 Table 2 Rearranged land use. New land use

87

Table 3 Impervious coefficients for each land use type. Before

Urban or built-up Airports; continuous and discontinuous urban fabric; roads, area port areas; quarrying; mining extraction zones and construction sites. Forest land Forest cover of: cork, holm oak, chestnut, oak, eucalyptus and pine Agricultural land Arable land; permanent crops; orchards; meadows Barren land Bare rock, beaches, dunes, sands and soils without vegetation Wetlands Watercourses, pounds and reservoirs.

Donigian (2002) states the coefficient of determination (R2) values range for assessing HSPF model performance. For daily flows, R2 is considered to be very good when it is greater than 0.8; good when it is between 0.7 and 0.8; fair when it is between 0.6 and 0.7, and poor when it is less than 0.6. For monthly flows, R2 is considered to be very good when it is greater than 0.85; good when it is between 0.75 and 0.85; fair when it is between 0.65 and 0.75, and poor when it is less than 0.65. Model validation was performed after the final parameter values obtained through calibration. Validation was performed for a 4 year period of observed values (2003e2006) and goodness of fit between observed and simulated values was reassessed. Error statistics described above also were computed for the model validation output data. However, the study relies upon a relatively short duration of recorded data, specifically on streamflow observed data. This may result in difficulty to identify trends in hydrological data, even if identified it may be difficult to attribute it to climatic variability due human intervention on land use or irrigation methods.

Land use type

Impervious coefficient

Urban or builtup area Agricultural land Forest land Barren land Wetland

0.6 0.4 0.3 0.1 1.0

most common value of 0.7, while all other parameters follow a uniform distribution. To explore the uncertainty in this parameter space (Refsgaard et al., 2007; Jeremiah et al., 2012), the GLUE method was implemented with repeated model simulations using HSPF hydrologic parameter values that were randomly selected from a predetermined probability distribution (Monte Carlo simulations) (K. Beven and Binley, 1992). This method was used to derive the 97.5% and 2.5% uncertainty bounds (95% confidence interval) for HSPF predictions of monthly discharge and evaluate the relationship between model uncertainty bounds and the threshold chosen for the likelihood measure (Lindblom et al., 2007a,b; Freni et al., 2008, 2009; Dotto et al., 2011; Vezzaro and Mikkelsen, 2012). To ensure that the Monte Carlo Simulations (MCS) adequately approximate model output uncertainty, 1000 model runs were performed. To calculate the model uncertainty, quantiles of time series weighted by the NasheSutcliffe coefficient (E) (2.5% and 97.5%) for a 95% confidence interval of discharge were derived by sorting ascending the sample discharges values at each month for behavioral parameter sets (where E was above or equal to 0.5). If x is period in months (January through December) and there are m samples (m ¼ 1e1000) and n are data months (n ¼ 1e60) we will get:

2.4. Uncertainty and sensitivity analysis Uncertainty analysis was conducted using a Monte Carlo approach. Multiple model simulation runs were performed using the 5 years period, with values for selected model parameters randomly chosen from assigned probability distributions. Table 4 shows the parameters used in the uncertainty analysis, their units and distribution. Minimum and maximum values were obtained from BASINS Technical Note 6. The parameter, IRC, follows a triangular distribution since it has a

xð1; 1Þxð1; 2Þxð1; 3Þ.xð1; nÞ; xð2; 1Þxð2; 2Þxð23Þ.xð2; nÞ; . xðm; 1Þ  ðm; 2Þxðm; 3Þ.xðm; nÞ; Sort ascending by the NasheSutcliffe Efficiency coefficient (for behavioral samples only):

Fig. 2. Lena River Basin.

88

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

Table 4 Parameter distribution, designation and calibrated values following model calibration. Name

Units

Range

Calibrated value

Designation

Function of

LZSNa AGWRCa DEEPFR

mm day1 none

(76.2,203.2) (0.92,0.99) (0.0,0.2)

76.2 0.97 0.1

Soil climate Baseflow recession Geology, ground water recharge

BASETP

none

(0.0,0.05)

0.02

AGWETP

none

(0.0,0.05)

0

INTFWa UZSNa LZETPAa LZETPUa LZETPFa LZETPBa IRCa INFILTa

none mm none none none none day1 mm h1

(1,3) (2.54,25.4) (0.5,0.7) (0.2,0.7) (0.6,0.8) (0.1,0.4) (0.5,0.7,0.7) (0.254,25.4)

Lower zone storage nominal Basic ground water recession rate Fraction of ground water inflow which will enter deep ground water and be lost Fraction of potential evapotranspiration, which is fulfilled only as outflow exists Fraction of remaining potential evapotranspiration Interflow inflow parameter Upper zone storage nominal Lower zone evapotranspiration parameter Lower zone evapotranspiration parameter Lower zone evapotranspiration parameter Lower zone evapotranspiration parameter Interflow recession parameter Index to the infiltration capacity of the soil

a

2 20.32 0.5 0.2 0.6 0.1 0.4 3.81

Riparian vegetation Marsh/wetlands Soils, topography, land use Surface soil conditions, land use Vegetation type/density, root depth Vegetation type/density, root depth Vegetation type/density, root depth Vegetation type/density, root depth Soils, topography, land use Soils, land use

Hydrology parameters proposed by AquaTerra (2004).

xð5; 1Þxð5; 2Þxð5; 3Þ.xð5; nÞ; xð2; 1Þxð2; 2Þxð2; 3Þ.xð2; nÞ; . xð9; 1Þxð9; 2Þxð9; 3Þ.xð9; nÞ xðm  3; 1Þxðm  3; 2Þxðm  3; 3Þ.xðm  3; nÞ; where: x(5,1) < x(2,1) < . < x(9,1) < x(m  3,1), this means that trial 5 has the lower E value (0.5) and trial m  3 has the higher E value of all trials. The percentiles of all samples for the first month are: xð5; 1Þxð5; 2Þxð5; 3Þ.xð5; nÞ according with percentile : Pð5Þ ¼ Eð5Þ=EðallÞ; xð2; 1Þxð2; 2Þxð2; 3Þ.xð2; nÞ according with percentile : Pð2Þ ¼ Pð5Þ þ Eð2Þ=EðallÞ; . xð9; 1Þxð9; 2Þxð9; 3Þ.xð9; nÞ; xðm  3; 1Þxðm  3; 2Þxðm  3; 3Þ.xðm  3; nÞaccording with percentile : Pðm  3Þ ¼ Pð9Þ þ Eðm  3Þ=EðallÞ; Pm where: EðallÞ ¼ i¼1 Ei : The next step was to find the discharge values for each month whose percentiles are 2.5% and 97.5%. To identify parameter sets as acceptable (behavioral) or unacceptable (nonbehavioral), the user must specify a threshold for the likelihood measure. In this study we considered behavioral parameter sets when the NasheSutcliffe coefficient was equal or above 0.5, and the relationship between model uncertainty bounds and the threshold chosen for the likelihood measure was evaluated. A sensitivity analysis was performed to study the responses of the hydrological parameters of HSPF model for our study area. To perform a multi-parameter sensitivity analysis, GLUEWIN software from Joint Research Centre, JRC, (Ratto et al., 2001) developed in MatlabÔ environment was used. This software provides visual sensitivity analysis through scatter plots and cumulative distributions proposed by Spear and Hornberger (1980). The outputs (objective function values, NasheSutcliffe coefficient) are classified into two groups: the 50% smallest values and the 50% highest values. Plotting the cumulative distributions for the two classes (high and low values) it is possible to perform a sensitivity analysis for a given parameter. The further apart the distributions are the more important the parameter is to that specific model and therefore more sensitive. Alternatively, if the parameter values do not split, then the parameter is unimportant and less sensitive. Assessment of HSPF water balance parameter sensitivity was performed to evaluate the impact of specific parameters on the output variables of interest. Specific key water quantity parameters used are described in Table 4. Since the sensitivity of model results to parameters in a specific watershed depends on the combined impacts of climate and watershed conditions (AquaTerra, 2004), the studies reported in the literature used for comparison have an average annual precipitation range from 800 to 3837 mm and catchment areas of 16e 6965 km2.

3. Results and discussion The hydrologic parameter set values for the expert system calibration shown in Table 4 are within the range of those presented in literature (USEPA, 2000). Table 5 shows the HSPExp criteria achieved in that calibration process (bold values indicate a fair calibration achievement). All criteria parameters were met

except for evapotranspiration with 36.4% error. The total runoff was overestimated by 4.8%. As stated, this is a data limited watershed and ET was not available at any of the meteorological stations for Lis watershed. The evapotranspiration error presented here is associated with the difference between simulated ET and derived ET (from observed maximum and minimum temperature). Daily (Fig. 3) and monthly (Fig. 4) hydrographs from the HSPExp model were plotted with the respective observed discharge from the Lis gaging station to determine the relative errors and to compare the timings of the flood peaks. The peaks correlate very well for most of the entire period. According to the literature, the mathematical statistical criteria show good results for the daily and monthly simulation according to the coefficient of determination criteria and the NasheSutcliffe coefficient of efficiency. Table 6 shows the results achieved for the statistical methods previously discussed for both the calibrated and validated model. For validation purposes, the calibrated model was run for 2003e 2006 without changing any parameter values. Simulated daily and average monthly, were compared with respective observed data values, Figs. 5 and 6 respectively. Even though the NasheSutcliffe coefficient of efficiency and the coefficient of determination increased, showing a very good approach, the simulated flow is overestimated by approximately 13%. An increase of almost 10% due to calibrated model, still according to literature for monthly simulation, shows good results. Thirteen parameter sets were included in the Monte Carlo parameter simulations. Table 7 shows the statistical distribution of these parameters. More than 90 percent of the 1000 MCS model runs produced values of E greater than 0.5, 49.8% percent had an E

Table 5 Output error percentage for the Lis River watershed model calibration. Flow component (units)

Simulated

Observed

% Error

Criteria

Total Runoff (mm) Total of highest 10% flows (mm) Total of lowest 50% flows (mm) Evapotranspiration (mm) Total storm volume (mm) Average of storm peaks (cm) Summer flow volume (mm) Winter flow volume (mm) Baseflow recession rate

1770.1 848.6 99.3 2367.8 312.9 11.8 79.6 985.8 23.6

1689.6 881.9 109.9 3723.64 288.4 11.7 78.1 1047.5 23.4

4.8% 3.8% 9.6% L36.4% 8.4% 0.8% 1.9% 5.9% 1.1%

10% 15% 10% 10% 10% 15% 10% 10% 1%

Bold values correspond to fair parameter calibration achievement.

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

a)

Table 6 Criteria statistical values achieved.

0

50

20 40 40 60 Observed Simulated Precipitation

20

80 100

10

120

Precipitation (mm)

30

140 0 1/1/1985

1/1/1986

b)

1/1/1987

1/1/1988

12/31/1988

12/31/1989

Observed

2000 1600 1200 800 400 0

0

500

1000

1500

Simulated Fig. 3. Daily observed and simulated streamflows for station 15E/03: a) time series with precipitation and b) scatter plot (R2 ¼ 0.716).

a)

Model calibration (1985e1989) Daily Average monthly Model validation (2003e2006) Daily Average monthly

E

R2

0.70 0.84

0.72 0.84

4.8 4.5

0.84 0.86

0.86 0.96

12.9 12.8

Dv(%)

greater than 0.7 and 125 parameter sets were associated with an E equal or greater than 0.8. The results indicate that the HSPF platform is well suited for the simulation of river discharge in the Lis River watershed, achieving a 0.70 and 0.84 NasheSutcliffe efficiency coefficients for daily and monthly calibration respectively. The threshold value for the likelihood measure to classify the model parameter sets as behavioral was determined as equal or greater than 0.5. To illustrate how sensitive the simulation of monthly discharge in the Lis watershed is related to the choice of the threshold value when using GLUE method with E as an efficiency of the likelihood measure, the model uncertainty bounds for thresholds of 0.5 and 0.8 were calculated. The 95% confidence intervals of model uncertainty due to parameter uncertainty with thresholds of 0.5 and 0.8 are compared in Fig. 7. The confidence interval of the simulated discharge changes significantly depending on the chosen threshold value. The lower the threshold value the wider the confidence interval of model uncertainty. The percentage of monthly average flows that fall outside the 95% interval confidence associated for a threshold of 0.5 and 0.8 is 4.8% and 6.0% respectively. Although a 5 year calibration and 4 year validation might be considered a short duration to ensure a wide range of wet and dry periods as well as storm events of varying intensity and duration, the results give an indication how

16 Observed

14

a) 40

Simulated

12

10 20

8 6 4 2

30 30 Observed Simulated Precipitation

25 20

40 50 60

15

70 80

Dec-89

Jul-89

Dec-88

Jul-88

Jan-88

Jul-87

Jan-87

Jul-86

Jan-86

Jul-85

10

Jan-85

0

0

35

10 Streamflow (m3 s-1)

Streamflow (m3 s-1)

89

5

90

0 1/1/2003

b)

12/31/2004

12/31/2005

100 12/31/2006

b)

500

1600

400

Observed

Observed

1/1/2004

300 200

1200 800 400

100

0 0

100

200

300

400

500

600

Simulated Fig. 4. Monthly observed and simulated streamflows for station 15E/03: a) time series and b) scatter plot (R2 ¼ 0.842).

0 0

250

500

750

1000

1250

1500

Simulated Fig. 5. Daily observed and simulated streamflows for validated model at station 15E/ 03: a) time series with precipitation and b) scatter plot (R2 ¼ 0.860).

90

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

a) 350

a) 16

Observed

16

200

Simulated

300

400

Fig. 6. Monthly observed and simulated streamflows for the validated model at station 15E/03: a) time series with precipitation and b) scatter plot (R2 ¼ 0.956).

the model might perform for an independent period having similar conditions. Sensitivity analysis (Fig. 8) clearly indicate that changes in five specific parameters (LZSN, AGWRC, UZSN, DEEPFR and INFILT) result in significant changes in the hydrology output (i.e. the model is most sensitive to these parameters), since the 50% high values are set apart from the 50% low values for the output (NasheSutcliffe weight). This means that by changing those parameters within the acceptable boundaries, model performance changes significantly. Conversely, the five parameters, BASETP, AGWETP, INTFW, IRC and, LZETP for all land uses have little effect on the output of the hydrology simulation in the Lis River watershed.

Table 7 Statistical parameter distribution. Minimum calibrated value

Maximum calibrated value

Mean

Standard deviation

76.285 0.920 0.000 0.000 0.000 1.003 2.548 0.501 0.200 0.600 0.100 0.507 0.481

203.074 0.990 0.200 0.050 0.050 2.999 25.391 0.700 0.699 0.800 0.399 0.699 24.849

137.964 0.956 0.101 0.025 0.025 1.999 14.111 0.601 0.451 0.700 0.245 0.631 13.054

36.313 0.020 0.057 0.014 0.015 0.590 6.389 0.057 0.145 0.057 0.085 0.047 5.330

Jul-89

Jan-89

Jul-88

Jan-90

Jan-90

100

Jul-89

0

Jan-89

0

Jul-88

0

Jan-88

4

Jul-86

100

8

Jan-86

200

Observed

12

Jul-85

Streamflow (m3 s-1)

Observed

95% Confidence Interval (0.8)

Station: 15E/03

300

Jan-85

b)

Jan-88

b)

Jul-87

Jan-07

Jul-06

Jan-06

Jul-05

Jan-05

Jul-04

Jan-04

Jan-03

Jul-03

0

Jul-87

0

50

Jan-87

4

Jan-87

100

Jul-86

150

8

Jan-86

200

Observed

12

Jul-85

250

Jan-85

Simulated

LZSN AGWRC DEEPFR BASETP AGWETP INTFW UZSN LZETPA LZETPU LZETPF LZETPB IRC INFILT

95% Confidence Interval (0.5)

Station: 15E/03

Streamflow (m3 s-1)

Streamflow (m3 s-1)

300

Fig. 7. Comparison of 95% confidence interval of monthly streamflow due to parameter uncertainty with different thresholds: a) 0.5 and b) 0.8.

The sensitivity analysis results reported in the literature show similar frequency occurrences on ranking HSPF water balance parameters the most sensitive in HSPF modeling (Fontaine and Jacomino (1997), AlAbed and Whiteley (2002), Jacomino and Fields (1997), Chung and Lee (2009), Chou et al. (2007), Abdulla et al. (2009), Donigian and Love (2007) and Kourgialas et al. (2008)). Table 8 shows the most sensitive parameters for each study included in this research. By adding the parameters according to the frequency observed in the hydrology calibration process, we established a ranking system to classify the most important parameters to calibrate water balance in HSPF. According to Table 9, INFILT and LZSN are the most sensitive parameters followed by AGWRC and UZSN which shows that HSPF calibration for water balance is most sensitive to function of soil and land use parameters. This clearly indicates that HSPF hydrology modeling is most sensitive to precipitation patterns, soil characteristics (LZSN, INFILT and UZSN) and land use characteristics (INFILT and UZSN). Table 10 summarizes climate and geographic information, about the studies presented in Table 8, which can be helpful to explain the relationship between parameter(s) sensitivity and climate and geographic conditions. These studies range from high to low catchment areas and average annual precipitation as well as different Soil and Land Use composition. From the relationship analysis one concludes that LZSN and INFILT parameters are generally sensitive as they appear in eight out of nine of the research papers presented. Although more information is required regarding soil composition it seems INFILT is considered a sensitive parameter in soils rich with loam rather than only one type of soil (sand, clay or silt).

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

Fig. 8. HSPF parameters sensitivity analysis in the Lis River watershed.

91

92

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

Table 8 Most sensitive HSPF water balance parameters reported in literature. AGWETP Fontaine and Jacomino (1997) (Normal Flow) Fontaine and Jacomino (1997) (Flood Flow) AlAbed and Whiteley (2002) Chung and Lee (2009) Chou et al. (2007) Abdulla et al. (2009) Donigian and Love (2007) Kourgialas et al. (2008) This study

AGWRC

BASETP

DEEPFR

x x

x x x x x x

x x x

x

x x

x

Parameters

Occurrences

Ranking

INFILT LZSN AGWRC UZSN DEEPFR INTFW LZETP AGWETP IRC BASETP

8 8 5 5 4 3 2 1 1 0

A A B B C D E F F G

4. Conclusions The hydrological behavior of Lis River watershed was modeled. The model prediction compared reasonably well to the observed data on daily and monthly flow rates and flow duration, resulting in a NasheSutcliffe efficiency coefficient of 0.70 and 0.84 respectively. Model validation demonstrates a good representation of the observed data. This is the outcome of graphical and statistical comparisons and measures of the model performance for daily and monthly streamflow.

INTFW

x x x

IRC

LZETP

x x x

x x

x

Table 9 Number of occurrences and parameter ranking.

INFILT

LZSN

UZSN

x x x x x x x x

x x

x x x

Modeling uncertainty and parameter sensitivity for the HSPF model using the GLUE method was evaluated in the Lis River watershed, with the aim to identify parameters that demand the greatest attention in light of limited data and where limited financial resources are often more intensive. The application of GLUE based on the NasheSutcliffe coefficient led to a good prediction uncertainty regarding the coverage of measurements by the uncertainty bands. The results provided in this paper show that even though monitoring climate conditions and streamflow are important in model performance, attention should be devoted to soil and land use data collection (LZSN and INFILT) since these data have been demonstrated herein as one important factor when characterizing river hydrology. Decisions may be taken directly from primary data measurements, derived statistics or the results of many steps of modeling, but it is the collected data that build the support for these decisions. A dataset is clearly of great value as it is inevitably collected through a huge commitment of time and money. Developing a hydrologic foundation based on this data for river segments throughout a region, for a selected time period long enough to represent climate variability, will provide water managers more time allocation on other water management decisions and new less relevant data to be collected.

Table 10 Watershed characteristics of the studies researched for the parameter sensitivity analysis. Study Fontaine and Jacomino (1997)

AlAbed and Whiteley (2002)

Catchment area (km2) 16

6965

Average annual PCP (mm)

Soil

Land use

1372

Silty/very fine loam

800e950

Predominant loam and sand

80% Forest 10% Grass 10% Urban 78% Agriculture 19% Forest 3% Urban 43% Urban 40% Forest 13% Agriculture 88% Forest 4% tea garden e 67% Forest 12% Agriculture 10% Urban 9% Wetland 2% Barren 58% Pasture 29.4% Crops 8.5% Forest 2.8% Urban 46% Agriculture 24% Forest 20% Barren 10% Urban

Chung and Lee (2009)

287

1325

e

Chou et al. (2007)

303

3495e3837

e

4000 5058

150e400 1092

Clay and clay loams silty loam Sandy till, clay and silt

Kourgialas et al. (2008)

130

500e2000

Clay loam

This study

176

855

Silt loam or loam (46%), sandy clay loam (40%) and sand, loamy sand or sandy loam (14%)

Abdulla et al. (2009) Donigian and Love (2007)

A. Fonseca et al. / Environmental Modelling & Software 51 (2014) 84e93

Acknowledgments This work is supported by project PEst-C/EQB/LA0020/2011, financed by FEDER through COMPETE e Programa Operacional Factores de Competitividade and by FCT e Fundação para a Ciência e a Tecnologia. André R. Fonseca acknowledges his doctoral fellowship (SFRH/BD/69654/2010) by the Portuguese Foundation for Science and Technology (FCT) and he also acknowledges Idaho State University for sponsoring his 2011e2012 visiting scholar position. V. Vilar’s acknowledges Ciência 2008 Program. References 2000/60/EC, C. D., 2000. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 Establishing a Framework of Community Action in the Field of Water Policy. Abdulla, F., Eshtawi, T., Assaf, H., 2009. Assessment of the impact of potential climate change on the water balance of a semi-arid watershed. Water Res. Manag. 23 (10), 2051e2068. AlAbed, N., Whiteley, H., 2002. Calibration of the Hydrological Simulation Program Fortran (HSPF) model using automatic calibration and geographical information systems. Hydrol. Process. 16 (16), 3169e3188. Ames, D.P., Michaelis, C., Anselmo, A., Chen, L., Dunsford, H., 2008. MapWindow GIS. In: Shekhar, Sashi, Xiong, Hui (Eds.), Encyclopedia of GIS. Springer, New York, pp. 633e634. AquaTerra, 2004. BASIN/HSPF Training Handbook. US EPA. ASCE, 1993. Criteria of evaluation of watershed models. ASCE task committee on definition of criteria for evaluation of watershed models of the watershed management committee, irrigation and drainage division. J. Irrig. Drain. Eng. 119 (3), 429e442. Baloch, M.A., Tanik, D.A., 2001. Application of BASINS/HSPF in the Koycegiz-Dalyan Watershed in Turkey: a Developing Country Case Study in Watershed Modeling. In: River Basin Management VI, WIT Transactions on Ecology and the Environment. Bergman, M., Green, W., Donnangelo, L., 2002. Calibration of storm loads in the south Prong watershed, Florida, using BASINS/HSPF1. J. Am. Water Res. Assoc. 38 (5), 1423e1436. Beven, K.J., 2001. Rainfall-runoff Modelling. Wiley Online Library. Beven, K., Binley, A., 1992. The future of distributed models: model calibration and uncertainty prediction. Hydrobiol. Process. 6 (3), 279e298. Bicknell, B., Imhoff, J., Kittle Jr., J., Jobes, T., Donigian Jr., A., Johanson, R., 2001. Hydrological Simulation ProgrameFORTRAN: HSPF Version 12 User’s Manual. AQUA TERRA Consultants, Mountain View, California. Borah, D.K., Bera, M., 2004. Watershed-scale hydrologic and nonpoint-source pollution models: review of applications. Trans. ASAE 47 (3), 789e803. Carrubba, L., 2000. Hydrologic modeling at the watershed scale using NPSM1. J. Am. Water Res. Assoc. 36 (6), 1237e1246. Chou, W.S., Lee, T.C., Lin, J.Y., Yu, S.L., 2007. Phosphorus load reduction goals for Feitsui Reservoir watershed, Taiwan. Environ. Monit. Assess. 131 (1), 395e408. Christiaens, K., Feyen, J., 2001. Analysis of uncertainties associated with different methods to determine soil hydraulic properties and their propagation in the distributed hydrological MIKE SHE model. J. Hydrol. 246 (1), 63e81. Chung, E.S., Lee, K.S., 2009. Prioritization of water management for sustainability using hydrologic simulation model and multicriteria decision making techniques. J. Environ. Manag. 90 (3), 1502e1511. Crawford, N.H., Linsley, R.K., 1966. Digital Simulation in Hydrology ’Stanford Watershed Model 4. Doherty, J., Johnston, J.M., 2003. Methodologies for calibration and predictive analysis of a watershed model. J. Am. Water Res. Assoc. 39 (2), 251e265. Donigian Jr., A.S., 2002. Watershed model calibration and validation: the HSPF experience. In: WEF National TMDL Science and Policy 2002, November 13e16, 2002. Phoenix, AZ. WEF Specialty Conference Proceedings on CD-ROM. Donigian, A.S., Crawford, N.H., 1976. Modeling Nonpoint Pollution from the Land Surface. US Environmental Protection Agency, Office of Research and Development, Environmental Research Laboratory. Donigian, A.S., Davis, H.H., 1978. User’s Manual for Agricultural Runoff Management (ARM) Model (Available from: the National Technical Information Service, Springfield VA 22161 as PB-286 366, Price codes: A 08 in paper copy, A 01 in microfiche. Report). Donigian Jr., A.S., Love, J.T., 2007. The Housatonic River watershed model: model application and sensitivity/uncertainty analysis. In: 7th International IWA Symposium on System Analysis and Integrated Assessment in Water Management, May 7e9, 2007. Washington, DC. WATERMATEX Proceedings on CD-ROM. Donigian, A.S., Huber, W.C., Laboratory, E.R., Consultants, A.T., 1991. Modeling of Nonpoint Source Water Quality in Urban and Non-urban Areas. Environmental Research Laboratory, Office of Research and Development, US Environmental Protection Agency. Donigian Jr., A., Bicknell, B., Imhoff, J., Singh, V., 1995. Hydrological simulation Program-Fortran (HSPF). In: Computer Models of Watershed Hydrology, pp. 395e442.

93

Dotto, C.B.S., Kleidorfer, M., Deletic, A., Rauch, W., McCarthy, D.T., Fletcher, T.D., 2011. Performance and sensitivity analysis of stormwater models using a Bayesian approach and long-term high resolution data. Environ. Model. Softw. 26, 1225e1239. Dunker, J.J., Melching, C.S., 1998. Regional Rainfall-Runoff Relations for Simulations of Streamflow for Watersheds in Du Page County. Illinois. U.S. Geological Survey Water-Resources Investigations Report 98-40335. EPA, 2000. BASINS Technical Note 6-Estimating Hydrology and Hydraulic Parameters in HSPF (EPA-823-R00-012). Fontaine, T.A., Jacomino, V.M.F., 1997. Sensitivity analysis of simulated contaminated sediment transport. J. Am. Water Res. Assoc. 33 (2), 313e326. Freni, G., Mannina, G., Viviani, G., 2008. Uncertainty in urban stormwater quality modelling: the effect of acceptability threshold in the GLUE methodology. Water Res. 49, 2061e2072. Freni, G., Mannina, G., Viviani, G., 2009. Urban runoff modelling uncertainty: comparison among Bayesian and pseudo-Bayesian methods. Environ. Model. Softw. 24, 1100e1111. Gallagher, M., Doherty, J., 2006. Parameter estimation and uncertainty analysis for a watershed model. Environ. Model. Softw. 22, 1000e1020. Hope, A., Stein, A., McMichael, C., 2004. Uncertainty in Monthly River Discharge Predictions in a Semi-arid Shrubland Catchment. British Hydrological Society, pp. 284e290. Hydrocomp, I., 1977. Hydrocomp Water Quality Operations Manual. Hydrocomp Inc, Palo Alto, CA. Jacomino, V., Fields, D., 1997. A critical approach to the calibration of a watershed model. J. Am. Water Res. Assoc. 33 (1), 143e154. Jakeman, A.J., Letcher, R.A., Norton, J.P., 2006. Ten iterative steps in development and evaluation of environmental models. Environ. Model. Softw. 21, 602e614. Jeremiah, E., Sisson, S.A., Sharma, A., Marshall, L., 2012. Efficient hydrological model parameter optimization with Sequential Monte Carlo sampling. Environ. Model. Softw. 38, 283e295. Jin, X., Xu, C.Y., Zhang, Q., Singh, V., 2010. Parameter and modeling uncertainty simulated by GLUE and a formal Bayesian method for a conceptual hydrological model. J. Hydrol. 383 (3e4), 147e155. Johnson, M.S., Coon, W.F., Mehta, V.K., Steenhuis, T.S., Brooks, E.S., Boll, J., 2003. Application of two hydrologic models with different runoff mechanisms to a hillslope dominated watershed in the northeastern US: a comparison of HSPF and SMR. J. Hydrol. 284 (1), 57e76. Kourgialas, N., Karatzas, G., Nikolaidis, P., 2008. Simulation of the flow in the Koiliaris River basin (Greece) using a combination of GIS, the HSPF model and a karsticesnow melt model. In: iEMS 2008: International Congress on Environmental Modelling and Software. Li, L., Xia, J., Xu, C.Y., Singh, V., 2010. Evaluation of the subjective factors of the GLUE method and comparison with the formal Bayesian method in uncertainty assessment of hydrological models. J. Hydrol. 390 (3e4), 210e221. Lian, Y., Chan, I., Singh, J., Demissie, M., Knapp, V., Xie, H., 2007. Coupling of hydrologic and hydraulic models for the Illinois River Basin. J. Hydrol. 344 (3e4), 210e222. Lowe, S., Doscher, R., 2003. Modeling of urban watersheds using basins and HSPF. J. Environ. Hydrol. 11. Lumb, A.M., McCammon, R.B., Kittle, J.L., Division, G.S.W.R., 1994. Users Manual for an Expert System (HSPEXP) for Calibration of the Hydrological Simulation ProgrameFortran. US Geological Survey, Reston, VA. Mapoteca, 2008. http://web.letras.up.pt/mapoteca/ (accessed June 2012). Martinec, J., Rango, A., 1989. Merits of statistical criteria for the performance of hydrological models. J. Am. Water Res. Assoc. 25 (2), 421e432. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I e a discussion of principles. J. Hydrol. 10, 282e290. http://dx.doi.org/10.1016/ 0022-1694(70)90098-3. Ratto, M., Tarantola, S., Saltelli, A., 2001. Sensitivity analysis in model calibration: GSA-GLUE approach. Comput. Phys. Commun. 136 (3), 212e224. Refsgaard, J.C., van Sluijs, J.P., Hojberg, A.K., Vanrolleghem, P.A., 2007. Uncertainty in the environmental modelling process e a framework and guidance. Environ. Model. Softw. 22 (11), 1543e1556. Spear, R., Hornberger, G.,1980. Eutrophication in peel inleteII. Identification of critical uncertainties via generalized sensitivity analysis. Water Res. 14 (1), 43e49. Sudheer, K.P., Lakshmi, G., Chaubey, I., 2011. Application of a pseudo simulator to evaluate the sensitivity of parameters in complex watershed models. Environ. Model. Softw. 26 (2), 135e143. USEPA, 2000. BASINS Technical Note 6. Estimating Hydrology and Hydraulic Parameters for HSPF. Vezzaro, L., Mikkelsen, P.S., 2012. Application of global sensitivity analysis and uncertainty quantification in dynamic modelling of micropollutants in stormwater runoff. Environ. Model. Softw. 27-28, 40e51. Zhang, J., Ross, M., Trout, K., Zhou, D., 2009. Calibration of the HSPF model with a new coupled FTABLE generation method. Prog. Nat. Sci. 19 (12), 1747e1755.

Further reading USEPA, 1999. BASINS Technical Note 5. Using HSPF with BASINS/NSPM. Zhang, W., Montgomery, D.R., 1994. Digital elevation model grid size, landscape representation. Water Resour. Res. 30 (4), 1019e1028.