Automatic Calibration Tool for Hydrologic Simulation Program-FORTRAN Using a Shuffled Complex Evolution Algorithm

Hydrologic Simulation Program-Fortran (HSPF) model calibration is typically done manually due to the lack of an automated calibration tool as well as the difficulty of balancing objective functions to be considered. This paper discusses the development and demonstration of an automated calibration tool for HSPF (HSPF-SCE). HSPF-SCE was developed using the open source software “R”. The tool employs the Shuffled Complex Evolution optimization algorithm (SCE-UA) to produce a pool of qualified calibration parameter sets from which the modeler chooses a single set of calibrated parameters. Six calibration criteria specified in the Expert System for the Calibration of HSPF (HSPEXP) decision support tool were combined to develop a single, composite objective function for HSPF-SCE. The HSPF-SCE tool was demonstrated, and automated and manually calibrated model performance were compared using three Virginia watersheds, where HSPF models had been previously prepared for bacteria total daily maximum load (TMDL) development. The example applications demonstrate that HSPF-SCE can be an effective tool for calibrating HSPF.


Introduction
While some hydrologic model parameters are measurable, others are either difficult to measure or represent some system process in such a way that physically determining the parameter value is not possible.Often, those parameters that are not directly physically based are calibrated.Calibration is the process of adjusting selected model parameters to minimize the difference between the simulated and observed variables of interest [1,2].Parameter calibration is necessary when using spatially-lumped hydrologic models like the Hydrological Simulation Program-FORTRAN (HSPF) [3,4].Model calibration may be performed manually, or the processes can be automated using an optimization algorithm [5,6].Manual calibration can be laborious and time consuming.On the other hand, an automatic model parameter calibration has the potential to be quicker and less labor intensive [5,[7][8][9][10][11].
The HSPF model is widely used to simulate hydrological processes and water quality in order to better understand and address a variety of water quality issues such as total maximum daily load (TMDL) development.In routine HSPF applications, the model is typically manually calibrated with initial parameter estimates and thoughtful adjustments [4,5,12,13].With HSPF, manual calibration assistance is provided by decision support software, the Expert System for the Calibration of HSPF (HSPEXP), which has been developed to provide guidance for parameter adjustment [1].However, even when an expert system is used, the results of a manual calibration are still often dependent on the modeler's experience and expertise.Thus, use of software like HSPEXP does not ensure calibration consistency across all users [5,8,9,14].
Several researchers have tried to calibrate HSPF using the Parameter Estimation (PEST) software tool [5,[14][15][16].However, the Gauss-Levenberg-Marquardt (GLM) search algorithm employed in PEST is not necessarily capable of locating a global optimum solution, and its performance is dependent upon an initial parameter set specified by the user [16,17].Consequently, there have been few applications of PEST in the field of surface water modeling [5].
Recent studies have tried to calibrate HSPF using random, sampling-based heuristic algorithms.Iskra and Droste [14] found that the random multiple search method (RSM) and the Shuffled Complex Evolution method (SCE-UA) could find a parameter set providing better model performance statistics than with PEST employing the GLM algorithm.Sahoo et al. [4] calibrated the hydrologic components of HSPF using a generic algorithm (GA), but it has been suggested that the GA required greater computing resources and time for parameter calibration than SCE-UA, making running the model less efficient [18][19][20][21][22].
The SCE-UA algorithm developed by Duan et al. [23] has been extensively tested in many hydrologic modeling studies, and it is now regarded as one of the most robust and efficient algorithms for parameter calibration [14,[18][19][20][21][24][25][26][27][28][29][30][31].Despite this, the SCE-UA algorithm has not been widely used in HSPF applications presumably because there is no tool developed to link the two together.
Parameter calibration using a sampling-based method like SCE-UA can benefit directly from the recent advances in computing resources and techniques.Particularly, the use of parallel computing has become more popular in hydrologic modeling because of its proven capability and potential [32][33][34].Although there exists a variety of parallel computing options developed for saving computational time, most of them are too complicated for use in routine modeling practices.Some computing software provides built-in or add-in parallel computing functions that hydrologic modelers can easily adapt for their own uses.Of them, "R", is an open-source program language and computing environment that supports parallel computing [35].
Previous studies examining auto-calibration for hydrologic models showed that the auto-calibration method did not always lead to successful calibration in terms of solution robustness and computational efficiency due to the limitations of the algorithm used [4,5,14,36].In this research, we have linked the HSPF model with the SCE-UA algorithm in a parallel computing framework supported by R (HSPF-SCE) with the purpose of providing an alternative and efficient tool for automated parameter calibration of HSPF.The new tool/approach was used to calibrate the HSPF models developed for three watersheds in Virginia.Output from the manual and auto-calibrated models are compared to demonstrate the performance of the HSPF-SCE calibration tool/approach.This paper presents a detailed description of the newly developed HSPF-SCE tool and exhibits its capability with example applications.

Hydrological Simulation Program-FORTRAN (HSPF)
The HSPF model is a process-based, continuous, spatially lumped-parameter model that is capable of describing the movement of water and a variety of water quality constituents on pervious and impervious surfaces, in soil profiles, and within streams and well-mixed reservoirs [37,38].Hydrologic simulation in the model consists of three modules: impervious land (IMPLND), pervious land (PERLND), and reaches, i.e., streams, rivers, and reservoirs (RCHRES).The IMPLND module represents impervious surface areas and simulates only surface water components.The PERLND simulates hydrologic processes happening on pervious surface areas, including infiltration, evapotranspiration, surface detention, interflow, groundwater discharge to stream, and percolation to a deep aquifer.The RCHRES module simulates hydraulic behavior of channel flow using the kinematic wave assumption.Details about simulation mechanisms of the model can be found in Bicknell et al. [37].

Shuffled Complex Evolution Method (SCE-UA) Algorithm
Heuristic optimization methods that adapt sampling-based, random-search approaches can be useful when an objective function is discontinuous and/or derivative information cannot be obtained since they do not require continuity and differentiability of the objective function surface [23].Many studies have demonstrated that heuristic optimization methods can provide answers close to the global optimum of the solution space [18,20,21,[24][25][26][27]29].Of the available heuristic optimization methods, the SCE-UA algorithm developed by Duan et al. [23] combines the simplex direct search method with strengths of three evolution algorithms including controlled random search, competitive evolution, and complex shuffling.The SCE-UA algorithm has been widely used in hydrologic modeling because of its sampling efficiency, which is attributed to combining the strengths of multiple optimization algorithms [23,39].In this study, the SCE-UA optimization algorithm was adapted as a calibration method for the newly developed tool, HSPF-SCE, because of its proven efficiency and ability to find the global optimum.

The R Software Package
R is an open-source software programing language and software environment for statistical computing and graphics [35], which was developed and implemented using the General Public License (GPL) that facilitates its public access [40].The capabilities of R are extended through user-created packages that develop specialized libraries and techniques [41].R also provides useful parallel computing capabilities which a user can apply to intensive computational tasks [42].Two existing packages in R were adapted for the development of the HSPF-SCE.The Latin Hypercube Sampling (LHS) package [43] was used to improve the efficiency of random sampling of the SCE-UA optimization algorithm, and the Snowfall package [44] was employed to increase computational efficiency of parameter calibration by running the HSPF model with multiple parameter sets at the same time in HSPF-SCE.

The HSPF-SCE Framework
In HSPF-SCE, the SCE-UA optimization algorithm is fully coupled with HSPF using R (Figure 1).HSPF-SCE transfers a pre-specified number of parameter sets sampled by the SCE-UA algorithm to HSPF and then reports objective function values calculated using simulated HSPF output back to the SCE-UA algorithm.Each parameter set in the initial population (all parameter sets) includes values for the parameters that are being used to calibrate HSPF.Initial calibration parameter values are selected from predefined, uniform distributions using a LHS method.The uniform parameter distributions are bounded by values provided in the US EPA HSPF guidance document Technical Note 6 [45].For the parameter optimization, the population of parameter sets is partitioned into several sub-groups or complexes.As the calibration proceeds, each complex "evolves" independently according to the competitive complex evolution (CCE) algorithm [46].The evolved complexes are combined into the next parameter set population.Then that population is re-partitioned, or shuffled, into new complexes based on the order of objective function values of each parameter set.The evolution and shuffling procedure iterations continue until a pre-defined stopping criterion is met.A more detailed explanation of evolution and shuffling procedures can be found in Duan et al. [46].Once parameter set values are determined, each parameter set in a population is incorporated into HSPF by means of changing the corresponding parameter values in the HSPF User Control Interface (UCI) file.HSPF is then run using each parameter set in the population.When the model runs are completed, HSPF-SCE calculates the value of the objective function.Then, the calculated result is fed to the SCE-UA routine as a basis to search for the next parameter set.Plots and statistics for evaluating model performance are developed outside the model in post-process.
In the SCE-UA algorithm, size of the population (number of parameter sets in this case) is determined as a function of the number of parameters being calibrated (N) and the number of complexes p as defined by Duan et al. [46].As the number of complexes increases, the chance of locating parameter sets satisfying the HSPEXP criteria increases, while computational efficiency decreases.In this study, based on preliminary analysis of the relationship between the time required to locate optimum and population size, the number of HSPF parameters that will be calibrated was set to 10 (N = 10), and the number of complexes p was set to 24 (p = 24).This yielded a calibration population size of 504 (population = p × (2N + 1)).The HSPF-SCE application developed here allows a user to change the criteria to stop the SCE-UA optimization iterations.In the application presented here, HSPF-SCE stopped searching parameters when the difference between the average of the lowest ten objective function values and the lowest objective function value returned for any given population of parameter sets was ≤1.5%.It should be noted that a discussion about how one might choose the most appropriate convergence criterion or the number of complexes for the purpose of improving the efficiency of the optimization process goes beyond the scope of this study.
As mentioned earlier, HSPF-SCE provides a parallel computing option when multiple processors (or cores and threads) are available.In this study, a four-processor Intel(R) Core(TM) i7 CPU 870@2.93GHzchip was used allowing for parallel computing and parameter calibration.When using the HSPF-SCE tool, the parameter sampling and data flow happen in R. The HSPF code is not altered when implementing HSPF-SCE.

Objective Function
When calibrating hydrologic models, the calibration objective function(s) are typically goodness-of-fit measures (e.g., coefficient of determination (R 2 ), Nash-Sutcliffe Efficiency (NSE)), with each assessing the degree of agreement between observed and simulated variables.Objective functions and the model performance criterion used to evaluate model calibration should be selected considering the objectives of any given modeling effort and the characteristics of the candidate objective function(s).Many previous studies have shown that using a single objective function may lead to unrealistic calibration results [5,7,8,47].Using multiple objective functions, and thus multiple measures of goodness-of-fit, may allow one to consider different aspects of fit between simulated and observed variables [5].As previously discussed, HSPEXP is a decision-support system that aids users who manually calibrate HSPF by offering expert advice about which parameters to adjust and how.HSPEXP guidance suggests the use of multiple objective criteria when assessing the adequacy of an HSPF hydrology calibration (Table 1).
Table 1.HSPEXP model performance criteria for hydrologic calibration of HSPF (revised from Kim et al. [5]).Kim et al. [5], using PEST, applied a single composite objective function that combined six sub-objective functions based on the HSPEXP calibration criteria.In this study, we adopted an objective function uniformly weighted with six performance measures so that multiple evaluation aspects could be considered simultaneously in a single objective optimization framework (Table 2).For the objective function formulation used in this study, the objective function value can range from 0% to 600%, with 0% being perfect agreement between the simulated and observed data.It should be noted that the purpose of this study was to develop and demonstrate a reliable and efficient tool for automatic calibration tool/approach, not developing and evaluating the most appropriate calibration objective function.

Table 2.
Objective functions used in the HSPF-SCE tool (revised from Kim et al. [5]).

Absolute error of daily flow
In Table 2, ( ) is the sub-objective function, θ is the parameter set, Θ is the feasible parameter range, Q is daily flow, EX is the fraction of time that stream flow equals or exceeds a specific flow rate, is the number of selected storm events, P is peak flow, is the number of summer and winter months, is the number of time steps in each month, is the number of time steps in each storm event, and is a weighting factor.

Study Watersheds
Considering the availability of existing HSPF models that had been manually calibrated for bacteria TMDL development [48][49][50], three watersheds located in the Ridge and Valley Physiographic Province of Virginia were selected for this study (Figure 2).The Piney River watershed drains 123 km 2 in Amherst County and Nelson County, and its predominant land cover is forest (79%), followed by pasture (10%), cropland (6%), and residential (4%).A National Weather Service Cooperative Weather station is located at the Montebello Fish Hatchery (COOP ID: 445690) within the watershed, and daily streamflow discharge has been measured at the watershed outlet (gauging station ID: 02027500) by USGS.Model calibration and validation periods were set to 1 January 1991 to 31 December 1995 and 1 January 1996 to 31 December 2000, in which 21 and 16 storm events were identified, respectively.The Pigg River watershed, which is mainly located in Franklin County, Virginia, drains 186 km 2 directly into Roanoke River.The dominant land use in the watershed is forest at 72%, followed by pasture (23%), cropland (3%) and residential (2%).The watershed has a NCDC Cooperative Weather Station (Rocky Mount, COOP ID: 447338) and a USGS gauging station (ID: 02058400) at its outlet.The hydrologic simulation of HSPF was calibrated and validated using streamflow measurements made at the USGS station between 1 September 1989 and 31 December 1995 and between 1 June 1984 and 31 August 1989, in which 29 and 23 storm events were identified, respectively.

Selection of Calibration Parameters
HSPF represents hydrologic and hydraulic features of a watershed using fixed and process-related parameters [51].Fixed parameters represent the hydraulic features of the drainage network and physical properties of the drainage basin, such as length, slope, width, depth and roughness of a watershed and areas covered by different soil types, land covers, and slopes.Process-related parameters are used to describe hillslope processes including rainfall interception, infiltration, runoff generation and routing, soil moisture storage, groundwater discharge into stream, and evapotranspiration [37,51].
Based on the HSPF model manual [45], sensitivity analysis [52], and the authors' professional experience, nine parameters were selected for calibration (Table 3).The value for one of the nine parameters, UZSN, was allowed to vary between the winter season and non-winter season.Applying different values of UZSN for winter and non-winter periods increased the number of calibration parameters (N) to ten.The same value for each calibration parameter was used for all PERLNDs with the exception of INFILT.Since INFILT varied by PERLND, INFILT was changed by a multiplier which retains differences between INFILT values.Table 3. Calibration parameters for hydrologic simulation of HSPF and their ranges [45].Possible ranges of parameter values found in the US EPA HSPF guidance document Technical Note 6 [45] were used to define the parameter space.For each of the three watershed models used here, the parameters not selected for calibration were fixed and left unchanged from the values that were used in the manually calibrated TMDL models.

Model Performance Evaluation
There is no firm consensus when it comes to acceptable hydrologic model performance measures; there is no one statistic that can be used to assess all aspects of model performance [38].Thus, it is often recommended that one use multiple performance statistics in conjunction with graphical/visual assessments and other qualitative comparisons rather than relying on a single quantitative metric [38,53].Having said that, most decision makers want definitive calibration targets or tolerance ranges [38].Several studies have proposed general target ranges for various metrics to evaluate model performance.Donigian et al. [54] provided HSPF model users with general guidance on model evaluation statistics, and Duda et al. [38] noted that the tolerance range of percent error should be considered so that the modeler and model-results consumer may make a more informed assessment of the model's performance.Moriasi et al. [53] suggested using performance statistics like NSE, percent bias (PBIAS) and RMSE-observations standard deviation ratio (RSR) and provided model evaluation guidelines for these measures (Table 4).A brief description of each Moriasi-suggested measure is provided below.Note: * Performance criteria ranges estimated from Figure 4 in Duda et al. [38].
R 2 describes the degree of collinearity between simulated and measured flow (Nagelkerke, 1991), ranging from 0 to 1, and is given by where N is the total number of flow data; Q is observed flow; Q is simulated flow; and the over bar denotes the mean for the entire evaluation time period.R 2 of 1 means a perfect linear relationship between two variables, while an R 2 of zero represents no linear relationship.
NSE is a normalized value that assesses the relative magnitude of the residual variance, ranging from minus infinity to 1 [55].NSE values greater than zero imply that the model predictions are more accurate than the average of the observed data, and a NSE = 1 indicates the model predictions completely match observed data.NSE is one of the most widely used statistics for assessing agreement between two variables in hydrologic modeling [53], and its use was recommended by ASCE [56] and Legates and McCabe [57].NSE is defined as PBIAS represents the overall agreement between two variables [58].A PBIAS of zero means there is no overall bias in the simulated output of interest compared to the observed data.Positive and negative PBIAS values indicate over-estimation and under-estimation bias of the model, respectively [58].PBIAS expressed as a percentage is given by Root mean square error (RMSE) is an absolute error measure commonly used in hydrologic modeling.Chu and Shirmohammadi [59] and Singh et al. [60] introduced RSR to facilitate relative comparison between RMSE values calculated for estimations in different units and scales by normalizing RMSE with the standard deviation of the observed data.RSR can vary from 0 to a large positive value, and a lower RSR value indicates better model performance [53].RSR is defined as In this study, calibrated HSPF hydrologic simulations were evaluated with statistical measures of R 2 , NSE, PBIAS, and RSR as wells as visual comparison of observed and simulated flow time series and flow duration curves.

Assessing Acceptable Estimated Parameter Sets
In this study, three HSPF hydrologic models were calibrated using the HSPF-SCE auto calibration tool.The HSPF-SCE tool was allowed to calibrate nine parameters, with one of those allowed to vary seasonally for a total of ten calibration parameters.In the calibration processes, the SCE-UA algorithm identified multiple parameter sets that satisfied the six HSPEXP model performance criteria while minimizing objective function values.For example, for the Reed Creek watershed, 252 parameter sets out of 504 possible parameter sets were found to meet all the six HSPXEP criteria.Figure 3 shows the Reed Creek distribution of the objective function values on the left y-axis and the number of HSPEXP criteria satisfied by the parameter sets in the last iteration of the optimization process on the right y-axis.Parameter sets began meeting all six HSPEXP criteria once the value of the objective function decreased to approximately 50%.The minimum objective function value was 11.6%.For the Piney River and Pigg River watersheds, 159 and 141 parameter sets satisfied all six HSPEXP criteria, respectively.The parameter sets that met all six HSPEXP calibration criteria are referred to herein as "qualified" parameter sets.Performance statistics produced by the qualified parameter sets are shown in Figure 4.For the Reed Creek watershed, hydrologic simulation using the 252 qualified parameter sets produced statistics in the "very good" ranges for monthly PBIAS, monthly NSE, and monthly RSR with some "good" measures for monthly R 2 (Table 4).The 159 qualified parameter sets for the Piney River watershed were in the ranges of between "good" and "fair" for monthly NSE and monthly RSR, "very good" for monthly PBIAS, but monthly R 2 values were in between "fair" and "poor".On the other hand, the 141 qualified parameter for the Pigg River watershed yielded relatively unsatisfactory performance statistics, and values of monthly NSE, monthly RSR, monthly R 2 and daily R 2 were classified as "poor".Karst topography including sink holes and springs that frequently appear in the Ridge and Valley physiographic region of Virginia [61] could be one possible reason for the poor model performance in the Pigg River watershed.HSPF has limited groundwater simulation capabilities, and representing karst hydrology using HSPF is challenging.
Once multiple parameter sets that met all the HSPEXP criteria were identified, a single parameter set expected to best represent the hydrologic processes of a study watershed was selected from the pool of qualified parameter sets.This selection was based on model performance statistics, visual comparisons of various model output graphics (e.g., Figure 5), and best professional judgment.For example, for the Reed Creek watershed, 160 parameter sets out of 252 were qualified, meaning they satisfied all six HSPEXP criteria for both the calibration and validation period.The 160 qualified parameter sets were then classified into groups based on five performance statistics as shown in Tables 5 and 6.Parameter sets belonging to the same group were regarded as equal in terms of model performance.In addition to the comparison of performance statistics presented in Tables 5 and 6, the qualified parameter sets classified into Group 1 were further assessed by visually comparing hydrographs and flow duration curve plots simulated using those Group 1 parameter sets.To illustrate the graphical/visual model performance evaluation, the daily flow time-series and flow duration curve simulated using two of the qualified parameter sets (i.e., No. 40 and No. 179 in Tables 5  and 6) are plotted in Figure 5.Although parameter set No. 179 provided better monthly RSR, PBIAS, and NSE than did No. 40 (Table 5), the graphical comparison (Figure 5) clearly shows parameter set No. 40 yielded a better match to the observed flow.Thus, parameter set No. 40 was selected as the final parameter set to simulate the hydrology for the Reed Creek watershed.The same parameter selection process was applied to the Pigg and Piney River watersheds.

Comparing Automated and Manual Calibration Parameter Sets
Manually calibrated parameter values were compared with the selected qualified parameter set identified by HSPF-SCE.Table 7 presents the comparison for all three watersheds, while Figure 6 illustrates those comparisons graphically for the Reed Creek watershed.The manual and automated approaches provided quite different ranges for some parameters: LZSN, UZSN, INFILT, BASETP, and AGWETP.Figure 6 shows box plots for the automated calibrated parameters for the Reed Creek model qualified parameter sets.In Figure 6, the interquartile ranges (IQR) of selected parameters (INFILT, AGWRC, DEEPFR, BASETP, AGWETP, IRC, INTFW, and UZSN-winter season) do not include the manually calibrated values implying that manual calibration is likely to fall in a local optimum in the parameter space.This finding does not agree with Kim et al. [5], who found general agreement among manually calibrated and PEST calibrated parameter values.The discrepancy between this study and Kim et al. [5] might be due to the use of different automated optimization algorithms (SCE-UA vs. PEST) and subjectivity in selecting a final parameter set from the pool of qualified parameter sets.

Comparison of Model Performance between Automatic and Manual Calibration
Hydrographs simulated with the selected qualified parameter sets were evaluated in terms of the six HSPEXP criteria.As seen in Table 8, both manually and HSPF-SCE calibrated parameters produced model output that meet all the criteria.The HSPF-SCE calibrated parameter set consistently provided lower bias in simulation of total volume compared to the manually calibrated parameters.Goodness-of-fit measures for the selected parameter sets are presented in Table 9.In general, the selected parameter values calibrated using HSPF-SCE provided performance statistics better than or equivalent to those calibrated manually.The measures of the Piney River watershed indicated "fair" to "very good" in the calibration period and "good" to "very good" in the validation period for both calibration methods.For the Reed Creek watershed, relatively great differences were found in the performance statistics compared to the other watersheds.HSPF-SCE provided statistics in the "very good" range, while those of the manual method were in "fair" to "good" in the both the calibration and validation periods.The selected parameter set for the Pigg River watershed gave "unsatisfactory" modeling results in terms of R 2 , RSR and NSE in the calibration period.In all the cases, PBIAS values fell in the "very good" range.
Piney River observed and simulated daily and monthly flow and flow exceedance curves are compared in Figures 7-9.Overall, flows simulated using the HSPF-SCE calibrated parameters are similar to those simulated using the manual calibration method, especially for baseflow and high peaks (Figure 7).In the Pigg River watershed, overestimation and underestimation of stream flow are found in 1992 of the calibration period and the first year of the validation period, respectively.The HSPF-SCE calibrated parameter values resulted in better simulation results under the low-flow conditions of 1986 and 1991 than the manually calibrated parameters.In general, the HSPF-SCE calibrated parameters   The time required to perform the calibration for the study watersheds was compared between the automated and manual calibration methods (Table 10).The computational time required by HSPF-SCE employing between one and four processors was documented for the Pigg River watershed.Parallel computing using two and four processors was 47% and 66% faster than using a single processor, respectively, which indicates that parallel processing is indeed more efficient.Comparing the time required for the automated and manual calibration, for the Pigg River watershed, manual calibration took 3.8 times as many hours when compared to the parallel processing time requirement.For the Piney and Reed Creek watersheds, manual calibration required 4.3 and 1.5 times longer than the automated calibration, respectively.The numbers of model runs required by HSPF-SCE and the manual method were relatively small for the Pigg River watershed and larger for the Reed Creek watershed, implying that parameter calibration was more difficult for the Reed Creek watershed than the Pigg River watershed.The manual calibration time spent was estimated based on data collected during the respective TMDL development projects.All calibrations were performed on an Intel 2.93 GHz quad core machine with 4 GB of RAM on Windows 8 in a 64-bit environment.The simulation time estimates shown in Table 10 for the HSPF-SCE account for computational time only.As presented here, there is an additional step that must be completed after HSPF-SCE has identified the pool of qualified parameter sets; this is the graphical comparison that must be performed by the modeler to select the final parameter set from the qualified parameter sets.The authors estimate that for each of the study watersheds presented here, this process of selecting the final parameter set from the qualified parameter sets took about one day (8 h).

Summary and Conclusions
An automated calibration tool for HSPF was developed, HSPF-SCE, and its capability/applicability was examined with existing HSPF models developed for three Virginia watersheds.Utilizing the R software environment, the new tool links the HSPF model to the SCE-UA optimization algorithm without any modification of the HSPF model.The R software environment also allows HSPF-SCE to utilize parallel computing resources, making the tool computationally efficient.HSPF models that had been previously assembled for bacteria TMDL development purposes in three watersheds in Virginia were calibrated using HSPF-SCE.Model performance for the auto-calibrated and manually-calibrated models was compared.
HSPF-SCE calibrated parameters outperformed the manually calibrated parameters in terms of model performance statistics and in terms of how long it took to calibrate the model (HSPF-SCE was quicker).HSPF-SCE identified multiple qualified hydrologic parameter sets satisfying all six HSPEXP criteria, suggesting HSPF-SCE can be an effective tool for hydrologic calibration of HSPF.Manually calibrated parameter values often fell outside of the IQRs developed using the qualified parameter set values, indicating the manual calibration method may fall in a local optimum in the parameter calibration space.It was also demonstrated that satisfying the HSPEXP criteria does not necessarily imply good model performance in terms of commonly used statistics such as NSE, R 2 , RSR, and PBIAS.
The applicability of the HSPF-SCE tool to efficiently and effectively calibrate the HSPF model was successfully demonstrated in this study.However, potential improvements remain.It is worth mentioning that since the tool itself could not recognize flaws in the HSPF model setup, e.g., erroneous FTABLEs, the model to be calibrated needs to be verified before using the HSPF-SCE tool to prevent "best fit" but improper modeling results.It should also be noted that selection of the most representative (final) parameter set from among the qualified ones relies on modeler experience and expertise.In addition, the optimization algorithm SCE-UA used in this study was developed for aggregated single objective function optimization, and there are times when multiple objective function aspects may need to be considered in hydrologic model assessment.For example, calibrating a model for bacteria TMDL development in Virginia requires a multi-objective optimization algorithm and framework.Although the aggregated single object function successfully identified multiple qualified parameter sets in the calibration, it could not provide the Pareto optimal surface, thus trade-offs between the sub-objective functions could not be examined.The continued development and testing of multi-objective function calibration for HSPF presents an interesting next step to study.

Figure 2 .
Figure 2. Locations of the study watersheds.The Reed Creek watershed is 703.5 km 2 in size and located in Wythe County of Virginia.The watershed mainly consists of forest (52%) and pasture/hay land (38%) with residential (8%) and cropland (2%).A National Climatic Data Center's (NCDC) Cooperative Weather station (Wytheville, COOP ID: 449301) is located 15 miles due west of the Reed Creek watershed outlet, where a USGS gauging station (ID: 03167000) is found.The hydrologic parameters of the HSPF model developed for the Reed Creek watershed were calibrated using streamflow measurements made between 1991 and

Figure 3 .
Figure 3. Objective function values and the number of HSPEXP criteria met by the final parameter set population using the SCE-UA algorithm of HSPF-SCE (Reed Creek watershed).

Figure 4 .
Figure 4. Qualified parameter set model performance plots.Data generated by running HSPF with each qualified parameter set, then comparing observed and simulated model output using four model performance measures (a) Monthly PBIAS and Monthly NSE; and (b) Monthly R 2 and Monthly RSR.The square, triangle, and diamond correspond to the parameter set selected by the authors for subsequent HSPF simulations and model performance evaluation.

Figure 5 .
Figure 5.Comparison of the observed and simulated daily hydrographs with parameter set No. 40 (a,c) and No. 179 (b,d) for the Reed Creek watershed; (a,b) are daily log-scale hydrographs; and (c,d) are flow duration curves of daily flow.

Figure 6 .
Figure 6.Comparison of parameter values calibrated by the HSPF-SCE and manually calibrated parameter values for the Reed Creek watershed (UZSN_W is for winter season and UZSN_N is for non-winter season).

Figure 7 .
Figure 7.Comparison of the observed and simulated daily and monthly hydrographs with the selected parameter set for the Piney River watershed.(a,b) daily hydrographs; (c,d) monthly hydrographs; and (e,f) flow duration curves of daily flow.

Figure 8 .
Figure 8.Comparison of the observed and simulated daily and monthly hydrographs with the selected parameter set for the Pigg River watershed.(a,b) daily hydrographs; (c,d) monthly hydrographs; and (e,f) flow duration curves of daily flow.

Figure 9 .
Figure 9.Comparison of the observed and simulated daily and monthly hydrographs with the selected parameter set for the Reed Creek watershed.(a,b) daily hydrographs; (c,d) monthly hydrographs; and (e,f) flow duration curves of daily flow.

Table 4 .
General guidance for performance assessment of hydrologic modeling.

Table 5 .
Classification of the parameter sets identified by HSPF-SCE in terms of model performance statistics (for the Reed Creek watershed).
Note: * A parameter set selected as the most representative at the final selection.

Table 6 .
The HSPEXP criteria values of Groups 1 to 4 (for the Reed Creek watershed).
Note: * A parameter set selected as the most representative at the final selection.

Table 7 .
Comparison of parameter values calibrated by HSPF-SCE and manually.

Table 10 .
Comparison of calibration time spent between automated and manual method.