Abstract

The U.S. Geological Survey developed the Expert System (HSPExp) in 1994 to assist modelers in the calibration of Hydrologic Simulation Program—FORTRAN (HSPF) models by calculating statistics, displaying graphs, and recommending parameter adjustments, based on model–data comparisons and several sets of rules. HSPExp has been used extensively since its release; however, HSPExp support was recently discontinued, and the software does not operate reliably on most modern computers. Many organizations have developed software to calculate calibration statistics, but these programs are neither freely available nor open source, which limits their usability by a wide user community. Using the existing open source codebase of Better Assessment Science Integrating Point and Nonpoint Sources (BASINS), RESPEC developed HSPEXP+ and released it in February 2015. HSPEXP+ replicates most of the functionality of HSPExp, and provides several additional reports and an advanced graphing functionality to assist a modeler in obtaining meaningful insights about the model on an overall basis, as well as on an operation-by-operation basis. This paper presents a case study of the calibration of the Snake River watershed situated in the eastern part of Minnesota. Snake River is a tributary to the St. Croix River at the border of Minnesota and Wisconsin, and drains a watershed of >2 500 km², with multiple lakes. HSPEXP+ was used to assist in hydrologic and water quality calibration in the watershed. The hydrologic model resulted in a model fit efficiency of >0.9 for monthly flow and >0.8 for daily flow with a mean error <10% in total volume. The water quality calibration included simulated loading rate comparisons with the target loading rates, and a graphical comparisons of instream observed and simulated data. The Snake River watershed model was satisfactory based on quantitative and graphical comparisons, and the model can be used for TMDL development.

1 Introduction

1.1 Background

The U.S. Geological Survey (USGS) developed HSPExp, the Expert System for calibration of HSPF (Hydrologic Simulation Program—FORTRAN) models in 1994 (Lumb et al. 1994). HSPExp was designed to assist HSPF users in calibrating the model by calculating calibration statistics, displaying graphs, and providing advice to improve the calibration by changing model parameters based on several sets of rules.

HSPExp interactively allows the user to edit the User Control Input (UCI) file for the HSPF (Bicknell et al. 2005) model, run HSPF simulations, produce plots of HSPF hydrology outputs as compared to observed values, compute error statistics for a simulation, and provide the user with expert advice on which parameters should be changed to improve the calibration. In general, use of HSPExp is an interactive process that requires repeating the cycle of simulation, computing statistics, reviewing plots, reviewing advice, and adjusting parameters.

HSPExp uses >35 hydrologic rules involving >80 conditions to recommend parameter adjustments to improve the flow calibration. The rules are divided into four groups: annual volumes, low flows, storm flows, and seasonal flows. Rules in a subsequent group are not tested until all rules in the previous group have been satisfied. The rules are based on the experience of experts in the use of HSPF in a range of climates and physiographic regions.

Output options include tables of statistics that compare observed data and simulated results, ten different types of graphical presentations, and tables of the expert advice. Time series output from the simulation is written to the watershed data management (WDM) file. Output may be viewed interactively or written to text files for printing or analysis in a spreadsheet program.

HSPExp has been used to support HSPF applications on numerous watersheds, such as the Connecticut watersheds (Love and Donigian 2002), the Housatonic River watershed (Donigian and Imhoff 2009) the LeSueur River watershed in Minnesota (Tetra Tech 2009). However, HSPExp was developed for older operating systems, and multiple users have reported that it is unsuitable for current computers. The USGS has stopped the support and active development of HSPExp.

1.2 HSPEXP+ Development

To replicate the statistics produced by HSPExp, many organizations developed in-house spreadsheets, scripts, and programs. However, these programs are neither widely available, nor are they open source or rigorously peer reviewed. Engineers of AQUA TERRA Consultants (acquired by RESPEC in September 2015) developed scripts to replicate HSPExp statistics and additional water balance calculations on the existing open source codebase of BASINS (Better Assessment Science Integrating Point and Nonpoint Sources, USEPA 2015a). The Department of Biological Systems Engineering (BSE) at Virginia Tech used these open-source scripts to develop the new interface and enhance its capabilities. RESPEC further improved the BSE program and released it for free with the name HSPEXP+ (Mishra et al. 2015) in February 2015. HSPEXP+ is open source software with the same licensing terms as BASINS. The entire codebase of BASINS and HSPEXP+ is available on Github, for users to view, clone, or fork. Since its release, RESPEC has been continuously updating and adding features to HSPEXP+ based on user input and client support. Although HSPF can simulate watershed models in traditional U.S. units or SI units, HSPEPX+ has been developed to support models developed with traditional U.S. units only, and therefore the HSPEXP+ outputs presented in this paper are in traditional U.S. units only.

RESPEC has used HSPEXP+ on multiple HSPF based watershed modeling studies (e.g. USEPA 2015b; Mishra et al. 2016). HSPEXP+ has also been downloaded at multiple universities and government organizations and is actively used for hydrology and water quality calibration. RESPEC is currently developing guidelines for the HSPF models developed for the Minnesota Pollution Control Agency (MPCA), so the models can be tested efficiently using HSPEXP+. This paper demonstrates the application of HSPEXP+ to the Snake River watershed model developed for MPCA.

1.3 Case Study

The U.S. Environmental Protection Agency (USEPA) requires the MPCA to carry out the Total Maximum Daily Load (TMDL) Program in Minnesota. In an effort to expedite the completion of TMDL projects, the MPCA has sponsored the construction of watershed models to support a simultaneous development of TMDL studies for multiple listings within a cataloguing unit or 8 digit Hydrologic Unit Code (HUC) watershed. In many cases, the MPCA is developing HSPF based watershed models. As part of the model development process, RESPEC was contracted to develop a watershed model for the Snake River watershed, shown in Figure 1 (HUC 07030004), along with other neighbouring watersheds (RESPEC 2016).

Figure 1 Location of Snake River watershed.

The Snake River watershed is mostly located in the Northern Central Hardwood Forest ecoregion of Minnesota. The watershed drains a land area >2 500 km², with forest, pasture, and wetlands as the major land categories. RESPEC used BASINS to develop the HSPF based watershed model for the Snake River watershed.

The meteorological data for the watershed was obtained using the BASINS network of stations and the local (Minnesota) volunteer citizen network. The latest concurrent meteorological data from BASINS was available from 1995 to 2009 for the stations covering the Snake River watershed. The atmospheric deposition data of ammonium and nitrate were obtained from the National Atmospheric Deposition Program (http://nadp.sws.uiuc.edu/). Only minor point sources existed in the watershed, and the point source data were provided by the MPCA.

The land use data for the watershed was obtained using the 2006 National Land Cover Dataset (NLCD), and the watershed delineation network was primarily obtained from the National Hydrography Dataset Plus (NHDPlus). The Snake River watershed had two long term USGS gauges: USGS 05335800 (Snake River near Pine City) and USGS 05337400 (Knife River). The USGS gauge on Snake River near Pine City was the primary calibration and validation gauge because it contained datasets for the complete simulation period 1995–2009, whereas the USGS gauge on the Knife River was discontinued after 2002. Observed water quality data as instantaneous samples were available at several locations in the watershed during the simulation period. The observed snow depth data were downloaded from the Midwestern Regional Climate Center’s cli-MATE website (http://mrcc.isws.illinois.edu/CLIMATE/) for several stations in and around the watershed. Based on the availability of data, the calibration period was set as 2002–2009, and the validation period was set as 1995–2001.

2 Results and Discussions

2.1 Snow Depth Calibration

Hydrology in Minnesota watersheds is significantly affected by snow because snow pack can remain for many months in the winter. Simulated and observed snow depths were plotted at several locations within the Snake River watershed for comparison. Observed snow depth data were available at several stations in and around the Snake River watershed. Although observed snow depth data were rarely available for the exact location that was being simulated, this situation for snow depth calibration is typical. These limitations do not allow a direct statistical comparison of observed and simulated snow depth, but a visual time series comparison can provide important insight into assessing the magnitude of snow depth, snow melt runoff, and timing of the snow melt period. These visual comparisons serve as a precursor to a complete hydrologic calibration.

HSPF can be set up to output time series of various modeled constituent concentrations and fluxes; however, plotting the individual time series using software packages like WDMUtil, GenScn, or BASINS4.1 is time consuming, and repetitive. HSPEXP+ allows the user to provide a graph specification file that specifies the output graph type (time series, frequency, scatter), location, axes details, data details, and display properties. For snow calibration, a time series of daily simulated average snow depth for forest and pasture land use categories was plotted along with observed data at the two nearest stations for each individual met segment (a region in the model that receives the same meteorological input) in the watershed (for example, see Figure 2).

Figure 2 Simulated and observed snow depth at met segment 50 (upstream watersheds) for years 2002–2004.

A snow depth frequency curve for the observed and simulated snow depth during the winter months was also prepared for each met segment (Figure 3). Visual analysis of these graphs for the 17 met segments suggests that the magnitude and timing of snow simulation were acceptable. No statistical comparison of observed and simulated snow depth was performed, as the observed snow depth data was not-continuous and the stations were several kilometres (5 km to 10 km) away from the met segments. Moreover, the hydrology (that implicitly includes snow simulation) calibration includes detailed statistical comparisons. Note that the hydrologic simulation is not strictly a linear process and, therefore, snow calibration was conducted in tandem with the hydrologic calibration.

Figure 3 Simulated and observed snow depth at met segment 50 for the calibration period.

2.2 Hydrologic Calibration

To conduct hydrologic calibration using HSPEXP+, a modeler needs to provide a BASINS specification file with an extension exs. This file provides HSPEXP+ with information about the location of the watershed, the WDM time series data file from which to access observed and simulated data, the data set ID numbers (DSN) of the observed and simulated data in the WDM file, storm dates, drainage area, seasons, and calibration criteria.

HSPEXP+ compares continuous observed and simulated daily flows at the calibration location and computes several expert statistics (adapted from Lumb at al., 1994) that include errors in total flow, flow percentiles, low flow recession, storm volumes, storm peaks, and seasonal flow. These statistics are output in a simple text file that users can refer to. Additionally, HSPEXP+ outputs error percent, Nash–Sutcliffe efficiency (NSE), root mean square error (RMSE), and RMSE-observations standard deviation ratio (RSR) (Moriasi et al. 2007) for total flow, yearly flow, monthly flow, and daily flow in simple text files.

Based on the expert statistics, HSPEXP+ generates advice on how to improve the hydrologic calibration (adapted from Lumb at al., 1994). HSPEXP+ follows a hierarchical order for analysis of statistics and advice generation based on whether the errors are within the limits of the criteria. The statistics are evaluated in this order: total volume, baseflow recession, high and low flows, storm peaks and volumes, and, finally, seasonal errors. The advice section is coded based on the advice provided in the original HSPExp program (Lumb et al. 1994). The user has flexibility in changing the error criteria, which allows skipping over certain statistics, if needed. Note that the advice section is limited, since it does not take into account snow or the variation of parameter values for different land uses.

Along with the reports that are generated in text files, HSPEXP+ also generates several time series graphs as portable network graphics (*.png) files that compare observed and simulated flow for the simulation period, and for the individual years, by using both arithmetic and log scales. These files are saved in date and time stamped folders for a user to refer back to. These graphs are in report quality format and can be directly imported to a model calibration report. Sample graphs and statistics generated by HSPEXP+ for the final calibrated Snake River model are shown in Figures 4 and 5 and Tables 1 and 2. The expert statistics suggested that the errors are within the criteria for most statistics, except 50% lowest flows and winter storm volume error. The winter volumes are generally difficult to calibrate in the Minnesota watersheds as the gauges may freeze in winter and the observed discharge in winter may be estimated rather than observed. The model fit efficiencies of > 0.8 for daily and >0.9 for monthly flows suggest a good to a very good calibration (Moriasi et. al. 2007).

Table 1 Part of expert statistics generated by HSPEXP+ for the Snake River model for the calibration period.

Statistics	Criteria	Statistic	Status
Error in total volume (%)	10	7.68	OK
Error in 10% highest flows (%)	15	–5.27	OK
Error in 25% highest flows (%)	10	–0.04	OK
Error in 50% highest flows (%)	10	6.25	OK
Error in 50% lowest flows (%)	10	20.09	Needs Work
Error in 25% lowest flows (%)	15	5.49	OK
Error in 10% lowest flows (%)	20	4.90	OK
Error in low-flow recession	0.03	0.01	OK
Error in storm volumes (%)	15	–1.43	OK
Seasonal volume error (%)	20	–10.82	OK
Error in average storm peak (%)	15	12.65	OK
Summer volume error (%)	20	–2.56	OK
Winter volume error (%)	15	8.26	OK
Summer storm volume error (%)	15	–7.54	OK
Winter storm volume error (%)	15	–18.62	Needs Work

Table 2 Hydrology statistics for the Snake River watershed at USGS gauge 05335800 for the calibration and validation period.

Monthly Flow Statistics	Calibration Period	Validation Period
Correlation Coefficient	0.91	0.94
Coefficient of Determination	0.82	0.89
Model Fit Efficiency	0.81	0.88
Daily Flow Statistics
Correlation Coefficient	0.87	0.86
Coefficient of Determination	0.76	0.74
Model Fit Efficiency	0.74	0.71

A complete hydrologic calibration is not only a statistical exercise in matching the observed and simulated data, but it also includes an analysis of the model water balance to verify that the model appropriately represents the properties of different land uses. HSPEXP+ generates water balance reports that present a detailed budget of the input and output of water in the watershed for each land use category as well as at every watershed outlet, and at annual and overall time intervals (for e.g. Table 3 below). These tables are useful in identifying the major pathways of runoff and the differences among the land use categories.

Visual and quantitative analysis of the calibration results suggested that the model satisfactorily represented the watershed characteristics and reproduced the flow at the calibration gauge. Following the hydrologic calibration, the model was run for the validation period. HSPEXP+ was used to produce hydrology statistics, graphs and water balance reports for the validation period. Some of the final validation statistics are presented in Table 2. After several iterations, the model was deemed acceptable for the calibration and validation periods, which meant it was ready for water quality calibration.

Table 3 Water balance (in.) report of Snake River watershed at the outlet for selected land uses.

Land Use	Forest - AB	Woody Wetlands	Cropland - AB	Watershed Total
Area (acres)	8,626	95,699	1,898	630,958
Influx
Rainfall	31.10	31.53	31.59	32.11
Runoff
Surface	0.12	0.01	0.44	0.50
Interflow	1.06	0.40	1.52	1.44
Base flow	7.31	7.76	7.97	7.80
Total	8.48	8.16	9.93	9.74
Groundwater Inflow
Deep	0.01	0.01	0.01	0.01
Active	7.78	9.89	8.05	8.62
Evaporation
Potential	33.22	32.37	34.57	33.52
Interception St	6.22	6.00	5.76	5.98
Upper Zone	4.79	4.73	5.33	5.20
Lower Zone	10.63	10.05	10.09	9.91
Groundwater	0.00	2.03	0.00	0.59
Base flow	0.37	0.14	0.13	0.22
Total	22.00	22.96	21.31	21.91

2.3 Sediment Calibration

Calibration of sediment and other water quality constituents starts by calibrating edge-of-stream sediment loading rates from each land use category. Similar to the water balance report, HSPEXP+ generates sediment loading reports (for example, Table 4 below). Observed data are rarely available to calibrate sediment loading rates; however, a literature review suggests target ranges of nutrient loading rates (AQUA TERRA Consultants 2015b) for different land uses. The sediment parameters were adjusted to ensure that the sediment loading rates are generally in the target rate ranges. The sediment loading rates calibration did not include statistical comparison as these loading rates are for general guidance, rather than the actual measurements in the watershed.

The instream sediment calibration process follows the loading rate calibration, and includes adjusting critical shear stress parameters for scour, and clay and silt deposition for each modeled stream reach (USEPA 2006). While transport of sand does not depend on the shear stress, its parameters should also be calibrated on a reach-by-reach basis. For the Snake River watershed, the computed shear stress in each reach was output to a WDM file at a 1 h time step. From these time series, 5th, 10th, 90th and 95th percentile exceedance values of shear stress for each reach were extracted. The 5th and 90th percentile values were used as initial values for critical shear stress for deposition and scour of clay, respectively and the 10th and 95th percentile values were used as initial values for silt. These values were later adjusted during the calibration process to ensure that the storages of sediment in the reaches remain generally stable.

Table 4 Sediment loading rates (tons/acre/year) for selected land uses in the Snake River watershed model.

Met Segment	Forest - AB	Forest - CD	Cropland - AB	Cropland - CD
Target Rates:	0.05–0.15		0.10–1.50
50	0.013	0.022	0.212	0.411
100	0.032	0.033	0.239	0.435
150	0.052	0.093	0.350	0.687
200	0.004	0.008	0.209	0.386
250	0.005	0.034	0.235	0.414
300	0.020	0.027	0.278	0.478
350	0.015	0.032	0.143	0.332
400	0.016	0.036	0.249	0.445
450	0.017	0.030	0.265	0.457
500	0.007	0.009	0.074	0.196
550	0.015	0.017	0.109	0.269
600	0.019	0.024	0.178	0.333
650	0.030	0.038	0.196	0.334
700	0.023	0.021	0.255	0.519
750	0.032	0.040	0.218	0.363
800	0.043	0.034	0.373	0.509
850	0.015	0.038	0.128	0.414
Mean	0.021	0.032	0.218	0.411
Max	0.052	0.093	0.373	0.687
Min	0.004	0.008	0.074	0.196

The instream calibration of sediment includes ensuring the general stability of each individual reach, unless site specific observation or data suggest otherwise. For example, excess scour may cause the bed depth to decrease rapidly leading to nearly complete depletion of the available storage of sediment for future storm events. It may also cause significant deposits of sediment in downstream water bodies. Unless, there is detailed river morphological data available to support a different behaviour, the reaches should be calibrated to ensure general stability. Again, this is not a statistical exercise so much as professional judgement.

The sediment budget report for each reach produced by HSPEXP+ was used to analyse the sources of sediment contribution to each reach and their general stability. The annual sediment reports for all reaches are more detailed, and show the deposition and scour of each size-fraction of sediment and the change in bed depth over time. The field visits and discussions with the local watershed coordinator suggested that the streams were generally stable in the Snake River watershed, and calibration was performed to maintain that stability. Table 5 below shows a portion of one of the sediment budget reports.

The stream sediment budget is observed in tandem with the instream concentration of sediment. The simulated instream concentration was plotted with the observed data (Figure 6) using HSPEXP+ by providing instructions in the graph specification file.

Table 5 Part of an annual sediment budget report for each operation.

Date	Mean	2002	2003	2004
Bed Depth (ft)	<NA>	5.02	5.021	5.021
Bed Storage (tons)
Sand	<NA>	112 570	112 690	112 780
Silt	<NA>	135 960	135 930	135 910
Clay	<NA>	47 248	47 203	47 172
Total	<NA>	295 770	295 820	295 870
Inflow (tons)
Sand	321.43	317.82	196.06	164.03
Silt	21 979	35 209	20 959	19 192
Clay	18 057	28 826	17 676	16 262
Total	40 357	64 353	38 831	35 618
Dep(+)/Scour(-) (tons)
Sand	213.11	199.76	122.6	99.345
Silt	−70.005	−70.909	−27.252	−17.994
Clay	−71.881	−66.114	−44.82	−30.806
Total	71.221	62.741	50.527	50.545
Outflow (tons)
Sand	108.32	118.05	73.46	64.683
Silt	22 049	35 280	20 987	19 210
Clay	18 129	28 892	17 721	16 293
Total	40 286	64 290	38 781	35 568

Figure 6 Observed and simulated total suspended solids at the outlet of Snake River watershed for the years 2002–2004.

Calibrating water quality to the observed data is extremely difficult, especially when there are no or very few data available during storm events, as in the case of the Snake River watershed. Observed total suspended solids (TSS) and other constituent data are available as a snapshot in time at a specific location and depth in a reach, whereas the model assumes the reach to be completely mixed; therefore the daily average data cannot exactly match this snapshot. However, comparisons with as much observed data as possible during different flow events helps in understanding the overall trends and ranges. Because of the paucity of data, quantitative statistical comparisons of instream observed and simulated TSS and other water quality constituents were not conducted. However, when HSPEXP+ plots the graphs, it also outputs the data used in the graphs to simple text files, which facilitates additional statistical analysis using other software.

In addition to the comparisons between observed data and simulated results, HSPEXP+ was also used to generate additional graphs, such as simulated TSS concentration at each reach and simulated detached sediment storage for all of the pervious land segments. Providing these instructions for hundreds of graphs in a graph specification file is simple, as these specification files can be produced in a spreadsheet program, and simple formulas can be used to generate instructions.

Although observed data is not available at every modeled reach, or for every land use category, these graphs help the modeler to quickly visualize the sediment behaviour in the entire watershed model. They also allow the modeler to quickly identify many types of problems in the model and ensure the reasonableness of the simulation. In the Snake River watershed, the visual analysis of graphs, visual and quantitative analysis of loading rates and reach budgets suggested that the model satisfactorily simulated TSS concentrations.

2.4 Nutrient Calibration

In the Snake River watershed, various forms of nitrogen (N) and phosphorus (P), their interactions and transformations, and other associated constituents—water temperature, dissolved oxygen (DO), biological oxygen demand (BOD) and phytoplankton—were modeled. The sources of these nutrients include point sources, nonpoint sources and atmospheric deposition.

The calibration of nutrients started with the calibration of loading rates of modeled nutrients, similar to sediment calibration. HSPEXP+ was used to generate reports of various forms of nutrients released from different parts of the watershed (Table 6 below), which gives the modeler an idea of the major pathways responsible for nutrient fate and transport in the waterways. Loading rates of BOD/organics, nitrate as nitrogen (NO₃-N), ammonia as nitrogen (NH₄-N), ortho-phosphorus as phosphorus (PO₄-P) were compared with the target loading rates developed for the midwestern United States (AQUA TERRA Consultants 2015b). During the loading rate calibration, an attempt is made to make sure that the loading rate is generally within the target ranges, but it is not a statistical exercise. After the simulated nutrient loading rates were adjusted to match the targets, the instream calibration process was started. Note that the loading rates calibration was revisited during the instream process calibration.

HSPEXP+ was set up to generate graphs of observed and simulated nutrient concentrations (for example, Figure 7 below). As with sediment, no statistical comparison of observed and simulated nutrient concentrations was conducted. However, visual analysis of graphs was used to verify that the simulated concentration was in the general range of the observed data. Instream nutrient calibration is generally conducted for all the nutrients in tandem, since most of the biochemical processes affecting these constituents are interrelated. For example, if the stream does not have enough biologically available N, phytoplankton stops growing and causes excess P to build up. Unless there is evidence supporting increasing or decreasing trends in nutrients or other related constituents, the model should be calibrated to ensure that the nutrient concentrations in the stream are relatively stable over time. HSPEXP+ was set up to generate graphs of multiple nutrients, phytoplankton as chlorophyll A, and benthic algae to visualize the interaction of nutrients in all the modeled streams and lakes in the watershed. This visualization allows the modeler to quickly ensure that all reaches in the watershed are behaving reasonably and to adjust parameters as needed (Figure 8 below). This graph is now available as a default option in HSPEXP+.

Table 6 Portion of the total nitrogen budget (lbs/acre/year) report at the outlet.

Land Use	Forest – AB	Woody Wetlands	Cropland - AB	Watershed Total
Area (acres)	8 626	95 699	1 898	630 958
NO3 (PQUAL)
Surface Flow	0.034	0.002	1.307	0.338
Interflow	0.026	0.009	0.99	0.228
Groundwater Flow	0.122	0.129	1.272	0.33
Total	0.182	0.14	3.569	0.896
NH3+NH4 (PQUAL)
Surface Flow	0.03	0.002	0.127	0.108
Interflow	0.017	0.006	0.093	0.039
Groundwater Flow	0.08	0.084	0.176	0.109
Total	0.127	0.092	0.396	0.256
Labile Organic N (PQUAL)
Surface Flow	0.014	0	0.124	0.054
Interflow	0.011	0.004	0.099	0.037
Groundwater Flow	0.063	0.067	0.913	0.16
Total	0.088	0.071	1.136	0.251
Refractory Organic N (PQUAL)
Surface Flow	0.031	0	0.281	0.121
Interflow	0.024	0.009	0.225	0.084
Groundwater Flow	0.142	0.151	2.07	0.363
Total	0.197	0.16	2.576	0.568
Total Nitrogen	0.594	0.463	7.677	1.971

In general, the simulated algae (phytoplankton and benthic algae) start growing in late spring, and the concentrations maximize in late summer. As algal growth increases, nutrient concentrations start to decline. Occasional spikes in nutrient concentrations are caused by runoff events that deliver nutrients from the land surface to the waterbody. The maximum growth of algae is limited by the growth parameters, availability of nutrients and other physical factors (light and temperature). The phytoplankton and benthic algae start to die off with colder temperatures and reduced light in the late fall. The nutrient concentrations start increasing as the algae die off. Deeper water bodies, like lakes, do not have significant benthic algae growth because of light extinction. If the model exhibits unusual behaviour, the modeler can go through the nutrient budget reports for the reaches to identify problems. The nutrient budget reports list all of the nutrients entering a waterbody by source as well as the instream gains and losses for each nutrient.

Figure 7 Observed and simulated nitrate as nitrogen plus nitrite as nitrogen at the outlet of Snake River watershed for years 2002–2004.

Figure 8 Concurrent plots of phytoplankton as chlorophyll a, benthic algae, nitrate as nitrogen, ammonia as nitrogen, and orthophosphate as phosphorus concentrations at the outlet of Snake River watershed for years 2002–2004.

HSPEXP+ also generates load allocation reports (AQUA TERRA Consultants 2015a) for sediment, TN and TP (e.g. Table 7). Following a satisfactory model calibration, these loading reports can be used to assist in total maximum daily load (TMDL) development and can help identify major sources of pollutants affecting the reach of interest. To calculate the effective contribution of each source in each reach, any losses of nutrients that may be happening in the reach are proportionally applied back to the sources. Any nutrient gain from the reach is listed as a separate source. This report illustrates the importance of each source based on its location in the watershed. A point source located upstream in a watershed will likely have a smaller effect on the nutrient loading at the outlet than a point source that is located closer to it.

Table 7 Load allocation report for total nitrogen in the Snake River and St. Croix-Kettle watersheds.

Source	Snake River Watershed
	Snake River–Cross Lake		Snake River at Outlet
	Load (lbs)	Percent	Load (lbs)	Percent
P:Forest–AB1	520	0.3	614	0.3
P:Forest–CD	19 842	10.4	19 998	10.1
P:Emerg Herb Wetland	2 608	1.4	2 638	1.3
P:Woody Wetlands	3 466	1.8	3 512	1.8
P:Grassland–AB	326	0.2	397	0.2
P:Grassland–CD	11 833	6.2	11 958	6.1
P:Pasture–AB	2 365	1.2	2 564	1.3
P:Pasture–CD	56 550	29.6	57 841	29.3
P:Cropland–AB	1 817	1.0	2 386	1.2
P:Cropland–CD	55 219	28.9	57 724	29.3
P:Cropland-Drained	15 552	8.1	15 996	8.1
P:Dev, Open Space	11 269	5.9	11 525	5.8
P:Dev, Low Intensity	1 969	1.0	1 973	1.0
P:Dev, Medium Intensity	601	0.3	600	0.3
I:Dev, Open Space	322	0.2	329	0.2
I:Dev, Low Intensity	271	0.1	271	0.1
I:Dev, Medium Intensity	377	0.2	376	0.2
Point Sources	5 570	2.9	5 920	3.0
Direct Atmospheric Deposition  on the Reach	0	0.0	0	0.0
Mass Balance Differences/  Additional Sources2	10	0.0	10	0.0
Diversion	0	0.0	0	0.0
Cumulative Instream Losses	–30 075	–15.7	–30 762	–15.6
Cumulative Instream Gains	608	0.3	608	0.3
Total3	191 096	100.0	197 243	100.0

¹P stands for pervious and I stands for impervious land uses.
²The additional sources may include sources other than nonpoint sources, point sources,atmospheric deposition, and upstream contribution.
³The total does not include losses because they have already been applied to the respective sources.

3 Conclusions

HSPEXP+ provides a modernized version of the original HSPExp program, plus many new functions that greatly facilitate both hydrologic and water quality calibration through generation of graphics, statistics, and complex summaries of HSPF output and observed data. This paper illustrates the application of HSPEXP+ to the Snake River watershed model that was developed using HSPF. The statistical and visual analysis of calibration results suggested that the hydrologic calibration and validation results were good to very good. For water quality calibration, no statistical tests were performed, but the quantitative analysis of loading rates and the visual analysis of observed and simulated concentrations suggest that the water quality calibration results were acceptable. Consequently, the Snake River watershed model can be used for TMDL development.

At the current stage, HSPEXP+ is an extremely useful and powerful tool to assist a modeler in conducting HSPF model calibration and validation. Application of HSPEXP+ can significantly reduce the resources required to generate useful peer reviewed statistics, reports, and graphs. As HSPEXP+ is open source and available for free, HSPF models developed by multiple agencies can be easily evaluated using same set of criteria.

Acknowledgments

The authors acknowledge the continued financial support of Clean Water Fund of Minnesota executed through the Minnesota Pollution Control Agency, St. Paul, Minnesota. The authors also acknowledge and thank the early adopters of HSPEXP+ for their feedback.

References

1. AQUA TERRA Consultants. 2015a. Documentation of Load Allocation Tool Developed under HSPEXP+. St. Paul, MN: Minnesota Pollution Control Agency. Project No. 21003-09.
2. AQUA TERRA Consultants. 2015b. Non-Point Source Target Loading Rates for Minnesota. St. Paul, MN: Minnesota Pollution Control Agency. Project No. 21003-09.
3. Bicknell, B. R., J. C. Imhoff, J. L. Kittle, Jr., T. H. Jobes and A. S. Donigian, Jr. 2005. HSPF Version 12.2 User’s Manual. Athens, GA:  U.S. Environmental Protection Agency, National Exposure Research Laboratory. Prepared by AQUA TERRA Consultants for Office of Surface Water, U.S. Geological Survey, Reston, VA and U.S. Environmental Protection Agency, Athens, GA. EPA/600/R-97/080.
4. Donigian, A. S., Jr. and J. C. Imhoff. 2009. Evaluation and performance assessment of watershed models. In Proceedings of the Water Environment Federation, TMDL 2009, 12–45. Alexandria, VA: WEF (Water Environment Federation).
5. Love, J., and A. S. Donigian, Jr. 2002. “A Modeling System of Nutrient Loads to Long Island Sound from Connecticut Watersheds.” Journal of Water Management Modeling 10: R208-11. doi: 10.14796/JWMM.R208-11.
6. Lumb, A. M., R. B. McCammon and J. L. Kittle, Jr. 1994. User’s Manual for an Expert System (HSPExp) for Calibration of the Hydrologic Simulation Program-Fortran. Reston, VA: U.S. Department of the Interior U.S. Geological Survey. Water Resources Investigation Report 94-4168. https://water.usgs.gov/software/HSPexp/code/doc/hspexp.pdf.
7. Mishra, A., A. S. Donigian, Jr. and B. R. Bicknell. 2016. HSPF Watershed modeling for the Upper St. Croix, Kettle, and Snake River Watersheds. St. Paul, MN: Minnesota Pollution Control Agency. Topical Report RSI-2604.
8. Mishra, A., P. B. Duda, M. H. Gray, B. R. Bicknell, A. S. Donigian, Jr. and R. W. Zeckoski. 2015. HSPEXP+ Version 1.11 User’s Manual. Mountain View, CA: AQUA TERRA Consultants.
9. Moriasi, D. N., J. G. Arnold, M. W. Van Liew, R. L. Binger, R. D. Harmel and T. L. Veith. 2007. “Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations.” Transactions of ASABE 50 (3): 885–900.
10. RESPEC. 2016. HSPF Watershed Modeling for the Upper St. Croix, Kettle, and Snake River Watersheds: Calibration and Validation of Hydrology, Sediment, and Water Quality Constituents. St. Paul, MN: Minnesota Pollution Control Agency. Report RSI-2606.
11. Tetra Tech. 2009. Minnesota River Basins Turbidity TMDL and Lake Pepin Excessive Nutrient TMDL. St. Paul, MN: Minnesota Pollution Control Agency. Report wq-iw7-32n. https://www.pca.state.mn.us/sites/default/files/wq-iw7-32p.pdf.
12. USEPA. 2006. EPA BASINS Technical Note 8: Sediment Parameter and Calibration Guidance for HSPF. Washington, DC: U.S. Environmental Protection Agency, Office of Water. https://www.epa.gov/sites/production/files/2015-08/documents/2006_02_02_basins_tecnote8.pdf
13. USEPA. 2015a. BASINS 4.1 Modeling Framework. Athens, GA: U.S. Environmental Protection Agency, National Exposure Research Laboratory, RTP, NC. http://www.epa.gov/exposure-assessment-models/basins.
14. USEPA. 2015b. Setup, Calibration, and Validation for Illinois River Watershed Nutrient Model and Tenkiller Ferry Lake EFDC Waters Quality Model. Washington, DC: U.S. Environmental Protection Agency. Submitted to USEPA Under Contract EP-C-12-052 Order No. 0002. https://www3.epa.gov/region6/water/npdes/illinoisriverwatershed/documents/final-watershed_lake_model_main_report%2008_07_15.pdf.