Data used to develop the HSPF model included: the 2006 land use data from the National Land Cover Database (NLCD) (
mrlc.gov/nlcd2006.php) [
12], elevation data from the National Elevation Dataset (NED) (
ned.usgs.gov/) [
13], and the hydrology and stream network from the USGS National Hydrography Dataset (NHD) (
nhd.usgs.gov/) [
14]. Soils data were obtained from the Soil Survey Geographic (SSURGO) [
15] database, including erosion potential (k-factor), texture (percent sand silt and clay), hydrologic soil group, available water storage capacity, and surface texture. Time series data including meteorological station data, streamflow, dam releases, snow water equivalent, and suspended sediment and sediment loads were collected for the time period 1958–2008. Station data for the meteorological development (air temperature, precipitation, and solar radiation) were collected from the National Weather Service (NWS) Cooperative Observer Program (COOP,
www.ncdc.noaa.gov/), Remote Automated Weather Stations (RAWS,
www.raws.dri.edu/), and California Irrigation Management Information System (CIMIS,
www.cimis.water.ca.gov/) stations within the study area. Solar radiation, wind speed, cloud cover and dew point were downloaded from the Fair Oaks CIMIS station. Streamflow, snow and sediment calibration and boundary condition data are described in the
Supplementary Materials (S1), the development of boundary dam conditions is described in
S2, and the development of hydrologic response units in
S3.
2.2.1. Gridded Meteorological Inputs to HSPF
HSPF does a reasonable job of predicting hydrology, sediment, water quality, and pollutant loads; however, the inability of meteorological stations to cover the spatial variations of precipitation and the uncertainty of managed flow data (diversions and return flows) are major limiting factors of more accurate predictions [
16]. A better spatial representation of watershed precipitation was found to improve model accuracy, even for a semi-lumped model [
17]. HSPF typically assigns data from meteorological stations directly to each sub-basin using Thiessen polygons. Meteorological stations are sparsely located in some areas and are not all active on any given day. Typical HSPF models rely on very few climate stations, some with only one station to represent the entire watershed [
18,
19,
20,
21]. The objective of the meteorological development was to enhance the spatial distribution of daily meteorological data and therefore increase the accuracy of modeled streamflow and suspended sediment. Incorporating gridded meteorological data as an input to HSPF has been successfully demonstrated [
17,
22] and, as the key driver of the rainfall runoff process and sediment transport, was shown to have increased the accuracy of the hydrologic results [
17].
Meteorological data from climate stations consisting of daily precipitation and air temperature were spatially interpolated over the model domain using the Gradient Inverse Distance Squared (GIDS) spatial interpolation [
23]. Spatial interpolation of meteorological data can be inaccurate when there are few stations and large distances between stations. The GIDS method provides accuracies at least as good as established kriging techniques without the complexity and subjectivity of kriging and the required station density [
23,
24]. For every active station, GIDS develops regressions for each day including the variables northing, easting, and elevation to interpolate to each grid cell, in this case, 270 m. This approach provides a detailed and localized incorporation of topographic and regional influences on precipitation (
Figure 2A).
Precipitation was summed and the temperature was averaged to transform the daily maps into monthly maps. PRISM data (Parameter-elevation Regressions on Independent Slopes Model) (prism.oregonstate.edu/) (
Figure 2B) are recognized as high-quality spatial climate data sets [
25] and were spatially downscaled using GIDS to the 270 m scale and used for comparison to match the measured data trend and to keep the regional monthly spatial structures intact. A ratio was developed for each grid cell using monthly PRISM values. The daily GIDS maps were multiplied by the PRISM ratio to produce daily meteorological values (
Figure 2C) that sum or average to exactly match PRISM monthly values. This method was applied to precipitation and air temperature data, and the data were spatially distributed to each HRU (
Figure 2D). This method is an improvement on the typical distribution of meteorological stations because the measured data spatial trend is preserved and the regional monthly structures are incorporated across the watershed to produce better interpolations when data stations are not present for any given day.
Potential evapotranspiration (PET) is an important component of the water balance equation in HSPF. PET can be estimated directly by HSPF or estimated using pre-processing methods and applied as a daily time series for each HRU and reach. For this project, estimates of daily PET were developed as a pre-processing step using the Priestley-Taylor equation [
26], simulated hourly solar radiation, including topographic shading and atmospheric parameters, estimated cloudiness, and the gridded daily air temperature estimates described above. Cloudiness was estimated using the Bristow-Campbell equation [
27] which was calibrated to CIMIS and RAWS daily measured solar radiation and the gridded air temperature data. This method of estimating PET is considered to be an improvement over internal HSPF calculations that use pan-evaporation measurements with an adjustment factor or other temperature-based equations to distribute PET to each HRU. The Priestley-Taylor equation is an energy balance equation that has a non-linear relationship of PET to temperature that is considered more appropriate when extrapolating PET into the future when elevated air temperatures are projected to persist [
28]. Other formulas that use a linear relationship of PET to air temperature have been shown to overestimate PET under warming climates [
28]. The PET estimates using the Priestley-Taylor equation were calibrated to all California Irrigation Management Information System (CIMIS) stations in the study area and then distributed as a daily time series to each HRU.
2.2.2. FTABLE Development
HSPF uses a hydraulic function table (FTABLE) to represent the geometric and hydraulic properties of stream reaches and reservoirs [
29]. FTABLEs define the stage-area-volume-discharge relation for each RCHRES. BASINS, the user interface that accompanies the HSPF code, generates FTABLEs from a simplified reach file using preset channel geometry assumptions. Some assumptions of the channel geometry include a trapezoidal-shaped channel cross-section, channel sides with slopes of a 1:1 ratio, the flood plain width on each side of the reach being equal to the mean channel width, and the maximum depth in the FTABLE being 100 times the depth of the floodplain slope break [
29].
One of the main problems with this simplified approach is that the BASINS-generated FTABLEs tend to have very few points within the range of the vast majority of actual flows. For example, in most of the HRUs there are only one to two data points in the BASINS-generated FTABLE that fall within all daily values for discharge. It is essential to develop FTABLEs that accurately reflect the major hydrological processes and water quality constituents such as sediment in order to ensure a realistic and unique simulation of the Sacramento River Basin. To enable a better simulation of suspended sediment, BASINS-generated FTABLEs were not used in the model; instead, FTABLEs were developed using hydraulic geometry. The heterogeneity of the Sacramento River watershed required two methodologies to develop FTABLEs; hydraulic geometry relationships were developed for the upland HRUs, and a one-dimensional hydraulic model was used for the lowland HRUs with very little slope.
Hydraulic Geometry Method
Regional hydraulic geometry relationships are often used to estimate channel dimensions for hydrologic models which require channel geometry data. Streamflow data from the National Water Information System (NWIS;
www.waterdata.usgs.gov) contains channel dimension measurements including the width, depth, and watershed size. These measurements were used to develop FTABLEs by ranking each gage by calculated discharge frequencies over a range of discharge exceedances. Analyses were limited to 45 out of the 288 gages that had greater than 10 years of data, had not moved or changed significantly during the period of record, and had less than 5% missing data. For each of the 45 gages, field measurements of channel width, depth, and velocity over a range of streamflows were converted to at-a-station hydraulic geometry parameters using a power-law fit for width, depth, and velocity with discharge:
where
w is the width (meters),
d is the depth (meters),
v is the velocity (meters per second),
Q is the discharge (cubic meters per second) and
a,
b,
c,
f,
k, and
m are the power-law fit parameters [
30].
The station power-law fit parameters (at-a-station) were scaled by each exceedance discharge (ranging from 0% to 100% of the maximum discharge). The power-law parameters for each exceedance discharge were then used to calculate the downstream hydraulic geometry [
30]. For example, at the 50% exceedance discharge for width, the at-a-station power-law fit was used to estimate the width for all of the gages that were used to calculate a single power-law fit for the downstream hydraulic geometry for the width at the 50% exceedance flow (R
2 = 0.97). Gages with zero discharge up to relatively high exceedances were removed because they violate the assumptions of the power law. Major reach segments (greater than a 10 km
2 drainage area and 10 km in length) were identified using GIS. Discharge was related to the drainage area for each gage and exceedance flow using a linear regression (R
2 = 0.93–0.97) (
Figure 3) because discharge was not known for the majority of reaches in the model domain.
The hydraulic geometry was calculated for each reach and exceedance discharge using the drainage area-to-discharge regression equations and the downstream hydraulic geometry relationship. Each HRU contains multiple reaches that were represented in the model as a single RCHRES; therefore, the reach segments were summed to generate the final FTABLE assigned to the RCHRES. Several reaches on the Sacramento River include flood structures called weirs. The weirs included in the HSPF model are located on reaches 61, 98, 99, 90 and 97 (
Figure 1). Adjustments were made to the FTABLEs to enable observed changes in streamflow routing during high flows. At certain high flow levels, the simulated weir will divert water and sediment to an adjacent bypass.
One-Dimensional Hydraulic Model Method
For seven of the lower-elevation HRUs in the valley floor (68, 78, 90, 94, 97, 98, 99), the hydraulic geometry method was not successful in formulating the FTABLEs. The breakdown of the hydraulic geometry method was primarily due to the low elevational gradients and large flood plains. For these lowland HRUs, the Digital Elevation Model (DEM)-derived drainage routing and watershed area frequently did not follow known river channels because out-of-levee floodplain areas are routinely near or below the river water level, and the steepest descent path did not accurately follow the river channels.
Instead of using the hydraulic geometry method for these low-elevation, low-gradient HRUs, an existing one-dimensional (1-D) hydraulic model was used to create the FTABLEs. The 1-D hydraulic model was created by the California Department of Water Resources (DWR) as part of the Central Valley Floodplain Evaluation and Delineation Program (CVFED) for the Sacramento River and its tributaries, with bathymetry and topography collected from 2001 to 2011 [
31]. The CVFED hydraulic model was run with flow only within the major tributary channels (no flows overtopped the levees). The hydrologic inputs to the model were the same as those used for generating the aforementioned hydraulic geometry FTABLES. The FTABLES were created from the 1-D hydraulic model by calculating the channel surface area, depth and volume for each reach over a range of streamflows.