Loading...
HomeMy WebLinkAboutFinal Neuse River Watershed Simulation Planrespec.com SIMULATION PLAN FOR THE NEUSE RIVER WATERSHED MODEL REPORT RSI-3348 PREPARED BY Cindie M. Kirby Paul Hummel Seth J. Kenner RESPEC 3824 Jet Drive Rapid City, South Dakota 57703 PREPARED FOR North Carolina Department of Environmental Quality Division of Water Resources 512 N. Salisbury Street Raleigh, North Carolina 27604 JULY 2023 Project Number W0392.22001.002 RSI-3348 i 2 TABLE OF CONTENTS 1.0 INTRODUCTION ............................................................................................................................................................... 1 1.1 PROJECT OVERVIEW ............................................................................................................................................................................ 1 1.2 PROBLEM STATEMENT AND BACKGROUND ................................................................................................................................. 2 2.0 DATA COLLECTION AND DEVELOPMENT ...................................................................................................................... 6 2.1 PRECIPITATION ...................................................................................................................................................................................... 6 2.2 EVAPOTRANSPIRATION AND OTHER METEOROLOGICAL DATA ............................................................................................. 8 2.3 STREAMFLOW DATA ............................................................................................................................................................................. 8 2.4 WATER QUALITY DATA ......................................................................................................................................................................... 10 2.5 POINT SOURCES .................................................................................................................................................................................... 12 2.6 ATMOSPHERIC DEPOSITION .............................................................................................................................................................. 14 2.7 OTHER DATA ........................................................................................................................................................................................... 15 2.7.1 Diversions and Withdrawals ................................................................................................................................................. 15 2.7.2 Irrigation.................................................................................................................................................................................... 17 3.0 SEGMENTATION AND CHARACTERIZATION .................................................................................................................. 18 3.1 DRAINAGE AREAS .................................................................................................................................................................................. 18 3.2 CHANNEL SEGMENTATION AND CHARACTERIZATION .............................................................................................................. 18 3.2.1 Reach Properties and Lake Selection ............................................................................................................................... 18 3.2.2 Numbering Scheme ............................................................................................................................................................... 19 3.2.3 F-Table Development ............................................................................................................................................................ 19 3.2.3.1 Lake F-Tables ............................................................................................................................................................ 19 3.2.3.2 Stream F-Tables ....................................................................................................................................................... 21 3.3 LAND SEGMENTATION ......................................................................................................................................................................... 21 3.3.1 Elevation ................................................................................................................................................................................... 22 3.3.2 Land Use ................................................................................................................................................................................... 22 3.3.2.1 Pervious and Impervious Land Classification ................................................................................................... 22 3.3.2.2 Septic Systems ......................................................................................................................................................... 24 4.0 CALIBRATION AND VALIDATION .................................................................................................................................... 26 4.1 CALIBRATION AND VALIDATION TIME PERIODS........................................................................................................................... 26 4.2 HYDROLOGY CALIBRATION AND VALIDATION PROCEDURES AND COMPARISONS ......................................................... 26 4.3 WATER QUALITY CALIBRATION ......................................................................................................................................................... 30 4.4 SENSITIVITY ANALYSIS ........................................................................................................................................................................ 32 5.0 DATA MANAGEMENT....................................................................................................................................................... 33 6.0 REFERENCES .................................................................................................................................................................... 34 RSI-3348 ii 2 LIST OF TABLES Table Page 1-1 Land-Cover Distribution ................................................................................................................................................................................ 4 2-1 U.S. Geological Streamflow Stations and Data Availability During the Model Calibration Time Period (January 1, 2004–December 31, 2022).................................................................................................................................................... 9 2-2 Pollutants That Will Be Calculated From the Point Sources ................................................................................................................. 13 2-3 Atmospheric Deposition Site Summary..................................................................................................................................................... 15 3-1 General Description of Hydrologic Soil Groups and Makeup of Watershed Areas ........................................................................ 23 4-1 General Calibration and Validation Targets or Tolerances for HSPF Applications ........................................................................ 29 4-2 Flow Calibration Criteria From Expert System for Calibration of HSPF ............................................................................................. 29 4-2 General Calibration Targets or Tolerances for HSPF Applications ..................................................................................................... 32 RSI-3348 iii 2 LIST OF FIGURES FIGURE Page 1-1 Neuse River Watershed Model Study Area ............................................................................................................................................... 3 2-1 Example Calibration Figure to Evaluate the Snowfall and Snow Depth Simulation ....................................................................... 7 2-2 U.S. Geological Survey Flow-Monitoring Station Locations ................................................................................................................. 11 2-3 Atmospheric Deposition Monitoring Station Locations ........................................................................................................................ 16 3-1 Land-Use Category Aggregation ................................................................................................................................................................ 23 4-1 Average Annual Precipitation for Modeled Watershed Areas from PRISM for Years 1990 to 2022 ......................................... 26 4-2 R and R2 Value Ranges for Model Performance ...................................................................................................................................... 28 RSI-3348 1 2 1.0 INTRODUCTION 1.1 PROJECT OVERVIEW The Simulation Plan for the Neuse River Watershed Model was developed to meet the requirements of SL 2020-18 Section 15(c) and calculate nutrient transport factors for the Neuse River Watershed with the watershed model and Scenario Application Manager. This watershed model will have a nitrogen- based focus and will also include phosphorus. The Neuse Nutrient Strategy (https://www.deq.nc.gov/about/divisions/water-resources/water- planning/nonpoint-source-planning/neuse-nutrient-strategy), implemented in 1997, addresses this problem by regulating major nutrient pollution sources throughout the Neuse River Watershed, including wastewater, urban stormwater, and agriculture. The strategy has succeeded partly by stemming additional nutrient loading during rapid population growth; however, the Neuse River Estuary remains impaired. The Neuse Nutrient Strategy is not supported by a calibrated watershed model, a standard tool for the development and evaluation of more modern nutrient strategies. A watershed model informs regulatory development and management decisions in many important ways, including by giving users a refined understanding of the influence of geography or various regulatory sectors on estuarine algal blooms. This project uses a contractor with oversight and support by the North Carolina Department of Environmental Quality (NC DEQ) Division of Water Resources (NC DWR) to develop a watershed model that meets agency standards for regulatory use and support in the Neuse River Watershed. The model will be a core product the NC DWR staff, stakeholders, and Environmental Management Commission rely on in their continual refinement of the Neuse Nutrient Strategy rules. This project will support long-overdue regulatory innovation that can drive systemic water quality improvements in the Neuse River Estuary. Recreation, property enhancement, recreational and commercial fishing, and greenways are some of the advantages of nutrient management. Excessive nutrient inputs can negatively influence the estuarine ecosystem and the communities that benefit from them. Conversely, these ecosystems and communities realize improvements from managing nutrient inputs. A watershed model is a critical scientific tool in structuring regulatory programs to achieve these broad-based environmental benefits, and several regulatory challenges or initiatives will be well -informed by the development of the Neuse River Watershed model. The Phase II Total Maximum Daily Load (TMDL), published in 2001 (https://www.deq.nc.gov/water-quality/planning/bpu/neuse/ neuse-tn-tmdl-ii/download), identified a specific need for watershed modeling. A dynamic watershed modeling approach is the most efficient means of obtaining detailed information on nonpoint-source and stormwater nutrient loads across the watershed. The dynamic watershed modeling approach will also help users better understand the impact of point sources. Directly measuring nutrient loads at the spatial and temporal scales of the watershed model would be impossible. Simplified watershed yield models provide annual nutrient loads; however, these models lack the temporal variability of loads, which is important for understanding episodic events or predicting loads under different climatic conditions. RSI-3348 2 2 The project team recommends the U.S. Environmental Protection Agency’s (EPA’s) Hydrological Simulation Program – FORTRAN (HSPF) [Bicknell et al., 2005 as the dynamic watershed model of choice for the Neuse River Watershed Model. Details of the model selection are presented in the model selection memorandum [Kenner, 2023]. Hydrologic Simulation Program-Fortran (HSPF) has been widely used throughout the United States to analyze water hydrology and quality in support of developing implementation plans based on attaining environmental goals. This complex and dynamic model can address soil, groundwater and surface-water processes, storm events, and impacts from point and nonpoint sources of pollution. The primary modeled parameters are nitrogen, phosphorus, suspended sediments, and streamflow. The model continues to be supported by the EPA and U.S. Geological Survey (USGS). The Neuse River Watershed Model study area is shown in Figure 1-1. 1.2 PROBLEM STATEMENT AND BACKGROUND RESPEC Company, LLC (RESPEC) proposes a project approach that will achieve the NC DWR’s project goals and objectives. The primary goal of this effort is to meet the requirements of Session Law (SL) 2020-18 Section 15(c) by developing nitrogen delivery factors for the Neuse River Watershed. The project will produce a calibrated watershed model for developing and evaluating modern nutrient strategies. As directed by SL 2020-18, the priority of this effort is to determine delivery factors for point-source discharges and nutrient offset credits. RESPEC will use the calibrated and validated watershed model to rigorously estimate the proportion of end-of-pipe nutrient loading from wastewater sources that reach the Neuse River Estuary. The delivery factors also play a key role in the availability and cost of nutrient trades between wastewater facilities and from wastewater facilities to watershed treatment best management practices; therefore, the watershed model we develop will provide a rigorous and unbiased approach to estimating relative nutrient contributions from the vast array of nonpoint nutrient sources throughout the watershed. The current Neuse Nutrient Strategy seeks to reduce nonpoint-source nutrients from cropland agriculture and new development while providing important protection by preserving riparian buffers. The watershed model can confirm existing nitrogen trading schemes and identify other nonpoint nutrient sources and their relative impacts on nutrient loading to the Neuse River Estuary. Understanding these sources offers opportunities to develop new trading strategies. The HSPF model is a comprehensive watershed model of hydrology and water quality that includes land-surface and subsurface hydrologic and water quality processes that are linked and closely integrated with corresponding stream and reservoir processes [Donigian et al., 2018]. HSPF is considered a premier, high-level model among those currently available for comprehensive watershed assessments and has experienced widespread usage and acceptance since its initial release in 1980, as demonstrated through hundreds of applications across the United States and abroad. HSPF is jointly supported and maintained by the EPA and USGS. HSPF is also the primary watershed model in the EPA BASINS modeling system and has been incorporated into the U.S. Army Corps of Engineers (USACE) Watershed Modeling System. This widespread usage and support have helped to ensure the continued code availability and maintenance for more than two decades despite varying federal priorities and budget restrictions. RSI-3348 3 2 Figure 1-1. Neuse River Watershed Model Study Area. RSI-3348 4 2 The Neuse River Watershed comprises four Hydrologic Unit Code 8 (HUC8) areas: the Upper, Middle, and Lower Neuse basins and the Contentnea Creek basin (Figure 1-1). These areas total approximately 6,062 square miles. The HSPF model being developed will characterize nutrient input from the Neuse River below Falls Lake, with a boundary condition established at the lake’s outlet using the USGS gage at that location (USGS 02087183, NEUSE RIV AT SR 2000 NR FALLS). The drainage area for this gage is 771 square miles, leaving a total area of 5,291 square miles in the model. Model results will be generated for this area; however, calibration will only be completed for the areas that are not tidally influenced. In coordination with NC DEQ, the model’s simulation period has been established as 2003 to 2022, with calibration beginning in 2004 after a warm-up period. Details of the model’s spatial and temporal scope are provided in the ensuing chapters. Table 1-1 shows the total area, average slope, and land cover, based on the 2019 National Land-Cover Database (NLCD) (https://www.mrlc.gov), of the project area (not including Falls Lake). The watershed has very diverse land use, with approximately one quarter crops, one quarter wetlands, one quarter forest, and one quarter other categories. Table 1-1. Land-Cover Distribution Watershed Neuse River Area (acres) 3,129,357 Average Slope (%) 2.7 Land Cover of Area (%) Open Water 1 Developed, Open Space 7 Developed, Low Intensity 5 Developed, Medium Intensity 2 Developed, High Intensity 1 Barren Land 0 Deciduous Forest 2 Evergreen Forest 14 Mixed Forest 6 Shrub/Scrub 2 Herbaceous 3 Hay/Pasture 3 Cultivated Crops 27 Woody Wetlands 24 Emergent Herbaceous Wetlands 3 The scope of work entails the following five major tasks: / Compile and Preprocess Data and Information to Support Model Development / Develop a Watershed Model of the Neuse River Watershed RSI-3348 5 2 / Apply Model to Establish Load Estimates / Develop Scenario Application Manager (SAM) / Deliver Model, SAM, and Documentation This report presents the simulation plan for developing hydrology and water quality model applications using HSPF for the Neuse River Watershed. This simulation plan presents the initial planned approach for constructing and calibrating the model applications with an emphasis on identifying and describing data requirements, sources, and availability. This plan will be revised after review comments are received, ongoing discussions with the participating agencies are completed, and the additional data needed to support the modeling effort are reviewed. This study plan, therefore, will be revised on an ongoing basis and will ultimately become part of the final report. The major steps in the model application development process consist of the following: 1. Collecting and developing time-series data 2. Characterizing and segmenting the watershed 3. Calibrating and validating the model These three steps are discussed in detail in the following sections. This report consists of seven chapters, including this introduction. Chapter 2.0 describes the collection and development of the hydrologic, meteorological, and other data needed for the simulations. Chapter 3.0 discusses the other types of spatial data needed for the segmentation and characterization of the watersheds. The calibration and validation process for the model is described in Chapter 4.0. Chapter 5.0 describes data management, and Chapter 6.0 includes references cited herein. RSI-3348 6 2 2.0 DATA COLLECTION AND DEVELOPMENT Hydrology and water quality simulation within the Neuse River Watershed requires the following types of time-series data: / Precipitation / Potential evapotranspiration / Other meteorological data (e.g., air temperature, wind, solar radiation, dewpoint, and cloud cover) / Streamflow / Water quality observations / Point sources / Atmospheric deposition / Other data (e.g., irrigation, diversions, and withdrawals) This section discusses the availability, selection, and processing methods of these time-series data for use in the watershed modeling. The quality assurance/quality control (QA/QC) and data management procedures are provided in Chapter 5.0. Only meteorological data are required to run the HSPF model; however, streamflow measurements and water quality observations are used to calibrate and validate the model. Other data types (e.g., point sources, atmospheric deposition, and diversions) help define the watershed’s inflow, outflow, and water quality. All the time-series data for the model will be placed into a Watershed Data Management (WDM) file, which is a binary database format originally developed to efficiently store large datasets used by HSPF and other models. 2.1 PRECIPITATION The HSPF model requires complete (i.e., no missing records) precipitation time-series data at an hourly timestep and with adequate spatial coverage and density across the model domains. Precipitation is the critical forcing function for all the watershed models because it drives the hydrologic cycle and provides the foundation for transport mechanisms that move pollutants from the land to the waterbody, where the pollutant impacts are imposed. The primary sources of long-term precipitation and other meteorological inputs for these watershed models include gridded data from the North American Land Data Assimilation System (NLDAS) and Parameter Elevation Regressions on Independent Slopes Model (PRISM). These data products are complete and available from 1979 to present (within the last few weeks of the download date). Because these data are gridded, they allow for easy extraction and aerial averaging over each hydrozone (i.e., an aggregation of subwatersheds that receive the same meteorological inputs) using scripted processes while providing efficient and consistent time-series extension. The NLDAS is a 12-kilometer (km) by 12-km dataset that provides hourly meteorological data. PRISM is a 4-km by 4-km dataset that provides daily precipitation totals, which are computed by combining a dense network of station data with radar measurement estimates that are interpolated based on a climate-elevation regression for each digital elevation model (DEM). Daily PRISM data will be used for RSI-3348 7 2 the modeling because these data provide a finer spatial resolution and generally a better fit to point precipitation data. The daily values will be disaggregated to an hourly timestep using the NLDAS data. The hourly NLDAS precipitation will also be loaded into the WDM file to provide another option to test during calibration. Specific stations are not associated with the gridded meteorological data. Precipitation data for the modeling period (January 2003 through December 2022) will be downloaded online from NLDAS (https://ldas.gsfc.nasa.gov/nldas/nldas-get-data) and PRISM (https://prism.oregonstate.edu/). The Neuse River Watershed typically gets less than 5 inches of snow per year. Because some snow falls occasionally, snow will be represented but will likely not be a significant portion of the calibration. Snow depth (i.e., snow on the ground) data are used to calibrate the snow accumulation and melt processes when the snow section of the model is active. These data are also used with mean and maximum winter- air temperatures to assess whether to activate the snow simulation capability within the watershed model. For the Neuse River Watershed, the snow depth (in inches) and snowfall (in inches) data are available through the National Climatic Data Center Global Historical Climatology Network stations [Menne et al., 2022] (ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily/). The snow depth data will be used during the hydrology calibration in multiple locations throughout the project area to ensure that snow processes are accurately represented. Graphs similar to that shown in Figure 2-1 will be developed for plotting snowfall, snow depth, and air temperature. Figure 2-1. Example Calibration Figure to Evaluate the Snowfall and Snow Depth Simulation. Snowfall and Depth for 398652 PERLND 552 Sno wfal l (in) RSI-3348 8 2 2.2 EVAPOTRANSPIRATION AND OTHER METEOROLOGICAL DATA In addition to precipitation, evaporation data are needed to drive the water-balance calculations in HSPF. Other meteorological time series are often required in temperate climates where snow accumulation and melt are significant components of the hydrologic cycle and water balance. These time series, such as air temperature (ATEM), solar radiation (SOLR), dewpoint temperature (DEWP), wind speed (WIND), and cloud cover (CLOU), are often required if soil and/or water temperatures are simulated. Water temperature is subsequently used to adjust rate coefficients in most water quality processes, and other time series are used in selected calculations (e.g., solar radiation affecting algal growth). The NLDAS dataset (https://ldas.gsfc.nasa.gov/nldas/) provides hourly ATEM, SOLR, and WIND parameters that will be directly applied to the meteorological time series with a conversion to the units needed for HSPF. The remaining meteorological constituents of CLOU, DEWP, and potential evapotranspiration (PEVT) are not directly available from the NLDAS dataset and require additional computations for this model. CLOU will be estimated for this model by SOLR data provided from the NLDAS database by using a parabolic equation [Thompson, 1976]. Two options for DEWP will be computed from a series of calculations that stem from the NLDAS specific humidity. The first option uses the specific humidity and ATEM to calculate the relative humidity [World Meteorological Organization, 2014]. Relative humidity will then be applied with ATEM to the August-Roche-Magnus approximation of the Clausius-Clapeyron equation [Stull, 2017] to calculate DEWP. The second option calculates a mixing ratio using specific humidity, and that mixing ratio is used with atmospheric pressure to estimate vapor pressure. DEWP is then calculated using the Clausius-Clapeyron equation [Stull, 2017]. Both options for DEWP will be assessed during calibration. Hourly PEVT estimates are included in the NLDAS dataset generated using a modified Penman [1948] energy-balance method; however, the NLDAS estimates of PEVT are included only for legacy compatibility with input requirements of the Sacramento Soil Moisture Accounting Model (http://hydromad.catchment.org/man/sacramento.html). The NLDAS PEVT estimates do not incorporate subsequent corrections to NLDAS estimates of energy forcing and have been found to overestimate evapotranspiration (ET) in other modeling efforts. Hourly PEVT will be represented by a computed Penman pan evaporation based on the Penman [1948] formula and the method of Kohler et al. [1955]. The necessary variables to compute the Penman pan evaporation are daily SOLR, DEWP, ATEM, and wind travel. Because two options for DEWP will be calculated, two options for PEVT will be calculated and assessed during calibration. 2.3 STREAMFLOW DATA Flow data are needed for calibrating and validating the watershed models to ensure that the hydrologic behavior of the watersheds, along with the transport of sediment and water quality constituents, is reproduced. Table 2-1 lists the continuous, observed streamflow data available in the Neuse River Watershed, with corresponding recording periods and the percent of missing data during the model calibration period (January 2004 through December 2022). The locations of the flow-monitoring sites RSI-3348 9 2 Table 2-1. U.S. Geological Streamflow Stations and Data Availability During the Model Calibration Time Period (January 1, 2004–December 31, 2022) (Page 1 of 2) Station I.D. Station Name Start Date End Date Missing (%) 02084909 SEVENMILE CREEK NR EFLAND, NC 06/24/1981 10/21/2004 95.7 02085000 ENO RIVER AT HILLSBOROUGH, NC 10/01/1927 01/30/2023 0 02085070 ENO RIVER NEAR DURHAM, NC 09/01/1963 01/30/2023 0.5 0208521324 LITTLE RIVER AT SR1461 NEAR ORANGE FACTORY, NC 09/30/1987 01/30/2023 0 0208524090 MOUNTAIN CREEK AT SR1617 NR BAHAMA, NC 10/01/1994 01/30/2023 0.1 0208524975 LITTLE R BL LITTLE R TRIB AT FAIRNTOSH, NC 10/24/1995 01/30/2023 0.1 02085500 FLAT RIVER AT BAHAMA, NC 08/01/1925 01/30/2023 0 02086500 FLAT RIVER AT DAM NEAR BAHAMA, NC 09/01/1927 01/30/2023 0 0208650112 FLAT RIVER TRIB NR WILLARDVILLE, NC 03/01/1988 09/29/2012 54 02086624 KNAP OF REEDS CREEK NEAR BUTNER, NC 10/01/1982 01/30/2023 10.7 0208675010 ELLERBE CREEK AT CLUB BOULEVARD AT DURHAM, NC 08/01/2008 01/30/2023 24.2 02086849 ELLERBE CREEK NEAR GORMAN, NC 10/01/1982 01/30/2023 10.7 02087183 NEUSE RIVER NEAR FALLS, NC 06/26/1970 01/30/2023 0 0208726005 CRABTREE CR AT EBENEZER CHURCH RD NR RALEIGH, NC 12/01/1987 01/30/2023 0 02087275 CRABTREE CREEK AT HWY 70 AT RALEIGH, NC 06/01/1997 01/30/2023 0 02087324 CRABTREE CREEK AT US 1 AT RALEIGH, NC 06/01/1990 01/30/2023 0 0208732534 PIGEON HOUSE CR AT CAMERON VILLAGE AT RALEIGH, NC 08/19/1987 01/30/2023 0 0208732885 MARSH CREEK NEAR NEW HOPE, NC 01/01/1984 01/30/2023 0 02087337 WALNUT CREEK AT BUCK JONES ROAD AT RALEIGH, NC 07/31/2018 01/31/2023 76.7 0208734210 WALNUT CREEK AT TRAILWOOD DRIVE AT RALEIGH, NC 08/08/2018 01/31/2023 76.8 0208734795 WALNUT CREEK AT SOUTH WILMINGTON ST AT RALEIGH NC 08/08/2018 01/31/2023 76.8 0208735012 ROCKY BRANCH BELOW PULLEN ROAD AT RALEIGH, NC 06/26/1992 01/30/2023 0 02087359 WALNUT CREEK AT SUNNYBROOK DRIVE NR RALEIGH, NC 05/01/1996 01/30/2023 0 0208739674 NEUSE R TRIB AT NRWWTP (CMP SITE) NR AUBURN, NC 05/09/2007 05/29/2008 94.4 0208739678 NEUSE R TRIB AT NRWWTP (CENTRAL STE) NR AUBURN, NC 06/07/2007 05/28/2008 94.9 0208741400 NEUSE R TRIB AT NRWWTP (EASTERN STE) NR AUBURN, NC 05/09/2007 05/28/2008 94.4 02087500 NEUSE RIVER NEAR CLAYTON, NC 08/01/1927 01/30/2023 0 02087580 SWIFT CREEK NEAR APEX, NC 03/01/2002 01/30/2023 0.3 0208758850 SWIFT CREEK NEAR MCCULLARS CROSSROADS, NC 12/01/1987 01/30/2023 0 0208762750 UNNAMED TRIB TO SWIFT CR NR YATES MILL POND, NC 05/01/2002 01/05/2011 89.6 0208762755 UNNM TRIB TO SWIFT CR AT NCSU RSRCH UNIT, RALEIGH 10/03/2009 01/05/2011 93.4 0208773375 SWIFT CREEK AT SR1555 NEAR CLAYTON, NC 10/01/2008 01/30/2023 25 02088000 MIDDLE CREEK NEAR CLAYTON, NC 01/01/2004 12/31/2022 0.1 02088383 LITTLE RIVER NEAR ZEBULON, NC 10/01/2008 12/31/2022 25.3 02088500 LITTLE RIVER NEAR PRINCETON, NC 01/01/2004 12/31/2022 0 RSI-3348 10 2 Table 2-1. U.S. Geological Survey Stations and Data Availability During the Model Calibration Time Period (January 1, 2004–December 31, 2022) (Page 2 of 2) Station I.D. Station Name Start Date End Date Missing (%) 02089000 NEUSE RIVER NEAR GOLDSBORO, NC 01/01/2004 12/31/2022 0 0208925200 BEAR CREEK AT MAYS STORE, NC 01/01/2004 11/29/2015 37.3 02089500 NEUSE RIVER AT KINSTON, NC 01/01/2004 12/31/2022 0 02090380 CONTENTNEA CREEK NEAR LUCAMA, NC 01/01/2004 12/31/2022 0 02090504 CONTENTNEA CREEK NR BLACK CREEK, NC 03/12/2021 12/31/2022 90.5 02091000 NAHUNTA SWAMP NEAR SHINE, NC 01/01/2004 12/31/2022 0 02091500 CONTENTNEA CREEK AT HOOKERTON, NC 01/01/2004 12/31/2022 0 0209173150 UNNAMED TRIB TO SANDY RUN AT SR1335 NR LIZZIE, NC 05/03/2006 02/17/2009 86.7 0209173190 UNNAMED TRIB TO SANDY RUN NEAR LIZZIE, NC 01/01/2004 02/18/2009 73.9 02091736 MIDDLE SWAMP NEAR FARMVILLE, NC 01/01/2004 03/14/2005 93.7 02091814 NEUSE RIVER NEAR FORT BARNWELL, NC 01/01/2004 10/24/2022 1 0209205053 SWIFT CREEK AT HWY 43 NR STREETS FERRY, NC 01/01/2004 06/30/2008 76.5 02092500 TRENT RIVER NEAR TRENTON, NC 01/01/2004 12/31/2022 0 are illustrated in Figure 2-2. As noted in Section 1.2, the Neuse River Watershed Model will be developed to characterize nutrient inputs for the area below Falls Lake, thus requiring a boundary condition to be established at the lake’s outlet. To quantify the Falls Lake boundary condition, streamflow and water quality data were downloaded from USGS gage 02087183 (NEUSE RIV AT SR 2000 NR FALLS). Flow data were downloaded from the USGS National Water Information System (NWIS) (https://waterdata.usgs.gov/nwis). All continuous streamflow data in the watershed will be included in the calibration; however, noncontinuous streamflow data are not as valuable for calibration purposes. As part of the model calibrations, the data will be plotted with a simulated flow. 2.4 WATER QUALITY DATA Water quality data are used primarily for model calibration and validation and to help quantify source contributions and boundary conditions. The specific constituents to be modeled in this study include all of the constituents needed for modeling nutrients. The following conventional constituents are modeled whenever nutrients are the purpose of a modeling effort: / Total suspended solids (TSS) / Water temperature / Dissolved oxygen (DO) / Carbonaceous biochemical oxygen demand ultimate (CBODu) (i.e., total CBOD) / Nitrite-nitrate (NO2/NO3) RSI-3348 11 2 Figure 2-2. U.S. Geological Survey Flow-Monitoring Station Locations. RSI-3348 12 2 / Total ammonia (NH3/NH4) / Total Kjeldahl Nitrogen (TKN) / Total nitrogen (TN) / Orthophosphate (PO4) / Total phosphorus (TP) / Phytoplankton as chlorophyll a / Benthic chlorophyll a Water quality data were collected from two sources: / USGS NWIS / NC DWR USGS NWIS water quality data were downloaded from the National Water Quality Monitoring Council Water Quality Portal (https://www.waterqualitydata.us/). NC DEQ provided data through 2021 from the North Carolina Ambient Monitoring System and Coalition Monitoring Program. These data will play a key role in the water quality calibration. For boundary condition, the Neuse River Near Falls, NC (USGS 02087183) gage has ample flow, and the NC DEQ water quality site Neuse Riv at SR 2000 nr Falls (J1890000) has ample water quality data to develop boundary condition loads with the FLUX software (https://water.usgs.gov/software/loadest/) for most parameters. Water quality parameters available at the J1890000 station include total ammonia, nitrate/nitrite, total Kjeldahl nitrogen, total phosphorus, total suspended solids, water temperature, and dissolved oxygen. Ratios of parameters throughout the watershed may be used for parameters not measured at the model boundary. Water quality data sources include the following: / USGS from Water Quality Portal (https://www.waterqualitydata.us/) / NC DEQ-provided water quality data 2.5 POINT SOURCES Point-source data throughout the Neuse River Watershed were provided by the North Carolina Department of Environment and Natural Resources (NC DWR). Discharging point sources with applicable data will be represented in the HSPF model applications. Applicable parameters available at the facilities include flow, 5-day biochemical oxygen demand (BOD5), ammonia nitrogen, nitrate and nitrite nitrogen, Kjeldahl nitrogen, total phosphorous, total organic phosphorus, total suspended sediment, dissolved oxygen, and temperature. The data, which were provided as daily data with values available approximately once a month (sometimes more or less), will be used to develop a filled daily time series following a set of rules and assumptions. Sites will be assumed as discharging every day of every month unless otherwise noted. Months with missing data can be filled in using the average of similar months (e.g., if January 2021 is missing data, the average of all of the other January data will be used to fill the month). Dates were RSI-3348 13 2 provided that represent when a permit was rescinded or operation ceased. In those instances, no flow will be represented after the date provided. It will be assumed that sites were operating before data were available because of the light nature of historical point-source data. Applicable parameters for the discharging facilities generally include CBOD5, ammonia nitrogen, Kjeldahl nitrogen, nitrate-nitrite nitrogen, DO, TP, TSS, and temperature. HSPF requires more input parameters than facilities typically provide or sample. A complete list of the required HSPF parameters is provided in Table 2-2. Table 2-2. Pollutants That Will Be Calculated From the Point Sources Pollutant Name Pollutant Description Daily Model Input Units Flow Effluent Flow Acre-Foot Heat Heat Energy of the Effluent BTU TSS Total Suspended Solids Tons DO Dissolved Oxygen Pounds NO3-N Nitrate as Nitrogen Pounds NO2-N Nitrite as Nitrogen Pounds NH4-N Total Ammonia as Nitrogen Pounds ORN Refractory Organic Nitrogen Pounds PO4-P Orthophosphorus as Phosphorus Pounds ORP Refractory Organic Phosphorus Pounds CBODu Ultimate Carbonaceous Organic Demand Pounds ORC Organic Carbon Pounds BTU = British thermal unit When facilities do not sample or report all the parameters listed, a dataset could be derived using a surrogate facility estimated with nutrient speciation factors or by setting a constant concentration, depending on the missing constituent. The assumptions that will be used for estimating missing parameters (provided in the following paragraph) have been applied to more than 50 HPSF model applications spanning several states and have been widely accepted by modelers, watershed managers, and point-source permitters. If dissolved oxygen or 5-day biochemical oxygen demand data are missing, concentrations of 8 and 1 milligrams per liter (mg/L) will be assumed, respectively. If nitrate-nitrite data are missing, a combined concentration of 2 mg/L will be assumed for non-wastewater facilities, and a combined concentration of 7 mg/L will be assumed for wastewater facilities. The combined nitrate-nitrite concentration will be partitioned into nitrate and nitrite based on other facilities with similar available data. If ammonia data are missing, a concentration of 2 mg/L will be assumed. Facilities without total phosphorus data will be assumed to have a total phosphorus concentration of 0.1 to 0.8 mg/L, depending on if the orthophosphate calculations result in a negative value; and 60 to 75 percent of the total phosphorus will be assumed as orthophosphate. Total phosphorus that is associated with biochemical oxygen demand RSI-3348 14 2 in HSPF is 0.7 percent, and the remainder of the total phosphorus that is not orthophosphate or associated with biochemical oxygen demand will be assumed as organic. Similarly, total nitrogen that is associated with biochemical oxygen demand in HSPF is 4.3 percent, and total nitrogen that is not nitrate, nitrite, or ammonia will be assumed as organic nitrogen. If total suspended solids data are missing, a concentration of 1 mg/L will be assumed, and all the total suspended solids concentrations will be split into 40 percent silt and 60 percent clay at each facility. Organic carbon will be assumed as 13 percent of the biochemical oxygen demand concentration. Besides temperature, concentrations of all the available constituents, including ultimate biochemical oxygen demand that will be converted from 5-day biochemical oxygen by using Equation 2-1 [Chapra, 1997], will be converted from mg/L to loads in pounds per day (lb/day) (i.e., concentration × flow × conversion factor; conversion factor = 8.34). Temperature will be converted from degrees Fahrenheit (°F) to a heat load in BTUs per day (i.e., temperature × flow × conversion factor; conversion factor = 8,339,145). −=−1 5 (5)1o k yLe (2-1) where: = = = u 55 1 CBOD CBOD 0.10(minimum value after primary treatment) oL y k Estimated daily time series will be imported into a WDM file, and loads will be applied to the corresponding stream in the external sources block of the user control input (UCI) file. 2.6 ATMOSPHERIC DEPOSITION Atmospheric deposition of nutrients is commonly included in watershed modeling efforts that focus on eutrophication issues. Nitrate and ammonium atmospheric depositions will be explicitly represented as a daily time series in the HSPF model applications. Wet atmospheric deposition data were downloaded from the National Atmospheric Deposition Program (NADP) (http://nadp.slh.wisc.edu/), and dry atmospheric deposition data will be downloaded from the EPA’s Clean Air Status and Trends Network (CASTNet) (https://www.epa.gov/castnet/). The nearest sites and corresponding recording periods are summarized in Table 2-3 and the locations are shown in Figure 2-3. Wet and dry atmospheric depositions will be applied evenly and directly to the waterbodies and land throughout the watersheds. The original dry deposition data are supplied at a weekly timestep as a particulate flux kilogram per hectare (kg/ha). The weekly data will be divided by seven to transform the data into a daily time series. The wet deposition is also supplied at a weekly timestep, but, in rare cases, sampling periods range from 1 to 8 days. Wet deposition data will not need to be divided by the number of days in the sampling period because wet deposition units are in concentration (mg/L) form. The concentration will instead be assigned to each day of the sampling period. In the model, the wet deposition data are multiplied by the precipitation amount to calculate the nutrient load. After being transformed to daily time-series data, the missing dry and wet deposition data will be filled in using interpolation when fewer than 14 missing RSI-3348 15 2 days have occurred between samples and by using monthly mean values when more than 14 missing days have occurred between values. The data will be converted to elemental concentrations and fluxes using multiplication factors from the UCI file (i.e., data are still nitrate and ammonia (not as N). The multiplication factors are used to convert the filled data into the units required by HSPF. The nitrogen deposition is applied as a time series to each segment, and the wash-off rates are mainly driven by precipitation intensity and calibration parameters. Table 2-3. Atmospheric Deposition Site Summary Site I.D. Name State Type Start Date End Date NC03 Lewiston NC Wet 10/31/1978 Present NC41 Finley Farm NC Wet 10/3/1978 Present NC35 Clinton Crops Research Station NC Wet 10/24/1978 Present NC29 Hofmann Forest NC Wet 7/2/2002 Present NC06 Beaufort NC Wet 1/26/1999 Present BFT142 Beaufort NC Dry 12/28/1993 Present CND125 Candor NC Dry 9/20/1990 Present Continuous wet and dry atmospheric phosphorus deposition data are not monitored through the NADP or CASTNet. Because of the lack of temporal data, an annual average value of total phosphorus deposition obtained from regional studies will be dispersed using the MONTH-DATA block in HSPF. Values of total phosphorus atmospheric deposition fluxes from the Chesapeake Bay area will be used for the model and range from 0.037 kilogram per hectare per year (kg/ha/yr) to 0.082 kg/ha/yr [Yang et al., 1996; Hu et al., 1998; Koelliker et al., 2004]. A midpoint value of 0.060 kg/ha/yr will initially be set, with higher values occurring in the summer and lower values occurring in the winter [Yang et al., 1996]. The total flux and monthly distribution may be adjusted as part of the calibration process. 2.7 OTHER DATA Additionally, ideal items to represent in the model application include groundwater and surface-water withdrawals, irrigation, and diversions information; these items would be represented using time-series data. If available, time-series data and/or estimates will be provided and processed to be included in the model. If time-series data are not available at a subwatershed or smaller level, estimations can be derived as described in the following sections. 2.7.1 Diversions and Withdrawals Time-series data for the surface-water and groundwater withdrawals were supplied by NC DWR. These values are available on a monthly basis and will be applied as a monthly value to represent stream diversion and withdrawals. Total surface-water withdrawals will be assigned to the reach they occur in and will be removed from the flow in each applicable reach. Individual withdrawals may be reduced as a part of the calibration process. Simulated withdrawals were also provided from the Oasis Model. Withdrawal data from NC DWR will be prioritized over Oasis simulated values. RSI-3348 16 2 Figure 2-3. Atmospheric Deposition Monitoring Station Locations. RSI-3348 17 2 2.7.2 Irrigation Irrigation in the Neuse River Watershed occurs on cropland, golf courses, and turf. Surface water and groundwater withdrawals for agricultural uses were supplied by NC DWR. As appropriate, these withdrawals can be applied to cropland in the hydrozone as precipitation if needed. RSI-3348 18 2 3.0 SEGMENTATION AND CHARACTERIZATION This section describes the methods proposed for the development of subwatershed, reach, and land-cover segments for the Neuse River HSPF model application. The segmentation and characterization define water travel from the various land uses within each subwatershed to each reach segment. 3.1 DRAINAGE AREAS Appropriate resolution for the subwatershed areas will be defined by the needs of North Carolina. Subwatersheds will be small enough to represent impaired reaches and lakes, as well as monitoring locations. The National Hydrography Dataset Plus Version 2 (NHDPlusV2) layers will be used to delineate watersheds where needed (https://www.nhdplus.com/NHDPlus/NHDPlusV2_data.php). Breaks will be made at HUC12 outlets, impairment outlets (impaired constituents), and monitoring locations. In most cases, significant point sources should be upstream of major monitoring and calibration locations within a subwatershed. The NHDPlusV2 dataset is a national, geospatial, surface-water framework that includes elevation, flow accumulation, and flow-direction grids. To delineate areas, batch points will be created in GIS at the desired breakpoints, and the ArcPro Watershed tools will be used with the NHDPlusV2. 3.2 CHANNEL SEGMENTATION AND CHARACTERIZATION The river channel network is the major pathway by which sediment and contaminants are transported from the watershed to each waterway; therefore, accurate representation or characterization of the channel system in the watershed for the model application is important. The river-reach segmentation considers river travel time, riverbed slope continuity, cross-section and morphologic changes, entry points of major tributaries, sampling locations, and impairment status. The channel characteristics are needed to define routing and stage-discharge behavior; bed composition for sediment, carbon, and nutrients; and bed/water-column interactions related to temperature, benthic oxygen demand, nutrient fluxes, and benthic algal mass. Because channel characteristics need to be defined spatially throughout the stream system, information from as many sites as possible will be used to define channel characteristics. 3.2.1 REACH PROPERTIES AND LAKE SELECTION The NHDPlus High Resolution flowline layer ((https://www.usgs.gov/core-science-systems/ngp/ national-hydrography/nhdplus-high-resolution) will be used to create the primary reach network. The primary reach layer was edited as needed by using the DEM and an imagery basemap. The North Carolina assessed streams and lakes (NC DWR Server – https://services2.arcgis.com/ kCu40SDxsCGcuUWO/ArcGIS/rest/services) will also be used as needed. The lakes selected to be explicitly modeled will be chosen based on impairment status, lake size, data availability, and location in the watershed. To be explicitly modeled, a lake was either impaired for a modeled parameter or was greater than 200 acres and connected to a modeled reach. RSI-3348 19 2 Reach length and slope are required to determine the physically based parameters in the model application and develop function tables (F-tables). These values will be calculated using ArcGIS for all non-lake reaches. Lakes that are modeled explicitly will be assumed to have an outflow; however, this assumption can be easily changed during calibration if any of the modeled lakes are determined as landlocked. Slope was derived from the USGS 10-meter (m) by 10-m, three-dimensional (3D) Elevation Program grid. 3.2.2 NUMBERING SCHEME This section describes the numbering scheme that will be used for the watershed drainage network. Reach identifications (I.D.s) are limited to one to three numerical digits in HSPF. Main-stem reaches occur along the main stem of each HUC8 watershed and will be given a reach I.D. that ends in zero (i.e., 0). These reaches will be assigned an odd-tens digit (i.e., middle number) if they represent a stream segment (e.g., 110, 130, 150, and 190 in the schematic) and an even-tens digit if they represent a lake (e.g., 120 and 160 in the schematic). Tributaries will be assigned an odd reach I.D. for the ones digit (i.e., end number) if they represent a reach (e.g., 141, 143, and 153 in the schematic) and an even number if they represent a reservoir (e.g., 142 in the schematic). The tens digit of the tributary reach I.D.s will correspond with the downstream, main-stem reach I.D. (e.g., 111 and 113 flow into 120). Reach I.D.s for subwatersheds and reaches will be numbered in order beginning with lower numbers upstream and ending with higher numbers downstream. If the next logical downstream, main-stem reach I.D. is not used, the downstream reach will be given the next largest main-stem reach I.D. For example, if a reach downstream of a main-stem reach with a reach I.D. of 170 and five tributary reaches (i.e., 171, 173, 175, 179, and 181) flow into the next downstream, main-stem reach, then that next main-stem reach would need to haven a reach I.D. of 190. Each subwatershed will typically only contain one waterbody (i.e., reach or lake) and will be given the corresponding reach I.D. 3.2.3 F-TABLE DEVELOPMENT This section describes the development of F-tables, which are required by the HSPF model to route water through each modeled reach (i.e., lake or stream). An F-table summarizes the hydraulic and geometric properties of a reach and is used to specify functional relationships among surface area, volume, and discharge at a given depth. 3.2.3.1 LAKE F-TABLES Data for lake F-table calculations include surface area and volume at various water elevations (depths) and overflow information. When available, surface area, volume, depth, spillway length, height above sill, and lake runout elevation data will be used for F-table development. Dam information was from the National Inventory of Dams (https://nid.sec.usace.army.mil/ords/f?p=105:1). Because these data are often unavailable, the F-tables will be based on the average values where data are missing, which is sufficient for the purposes of this model. If additional data become available, the data will be incorporated into the existing model application. The equations that will be used to calculate flows from lakes at different water elevations, as well as any assumptions made, are discussed in this section. For simplicity and because of the lack of overflow data, the equation of discharge for overflow spillways will be used to calculate discharge from lakes (Equation 3-1). Because of the project scale, coefficient correction factors for overflow calculations will not be used, and side contractions of the overflow and approach velocity have been disregarded, which allows for using the equation in its simplest form. RSI-3348 20 2 =1.5 eQ C L H (3-1) where: = Discharge (cubic feet per second[cfs]) = Water depth above weir (head, feet [ft]) = Effective length of crest (ft) = Variable coefficient of discharge. e Q H L C The total head ()H used in the equation will be calculated at variable water levels as the difference between the water-surface and outlet elevations. The outlet will be assumed at the maximum recorded depth (if available) or maximum contour depth. An effective length of the crest ()eL can be derived from a spillway length. When a spillway length is not available, the mean length of all of the available sites will be assumed. At lake depths below the outlet, eL will be set equal to the spillway length. At lake depths above the outlet, eL varies as a function of depth and will be increased assuming a 0.02 floodplain slope at each end of the crest. The variable coefficient of discharge ()C will be calculated using an empirical relationship derived by plotting x-y points along a basic-discharge coefficient curve for a vertical-faced section with atmospheric pressure on the crest from the U.S. Bureau of Reclamation [1987] (Equation 3-2): =  + 0.1528 3.8327 d PC In H (3-2) where: = Crest height (ft) = Head (ft). P H The crest height ()P will be assumed as the height above the sill (if available). The head ()H will vary with the water surface and will be calculated as described in the previous paragraph. When the height above the sill is unavailable, the mean value from all of the available sites will be assumed. After the available data are collected and combined, an F-table will be developed for each lake by calculating the surface area, volume, and discharge over a range of depths. F-tables for lakes will be developed using the calculated surface area, volume, and depth relations. For these lakes, the volume and surface area at incremental depths will be estimated using conical geometry and assuming a flat bottom for an inner circle with half of the radius of the maximum surface area. The highest contour (if available) or maximum depth will be assumed as the outlet. Depths will be added incrementally above the outlet until the F-table discharge exceeds the maximum observed discharge levels. The surface area and volume above the outlet will be calculated using conical geometry with an initial floodplain slope of 0.01. The discharge at each height above the outlet will be calculated using Equations 3-1 and 3-2. The discharge values of depths at or below the outlet will be zero. The initial value of the floodplain slope is arbitrary and can easily be adjusted during the calibration process. RSI-3348 21 2 3.2.3.2 STREAM F-TABLES Data requirements for stream F-table development include cross-section and discharge measurements. Cross-section measurements have been obtained from the North Carolina Flood Risk Information System HEC-RAS models (https://fris.nc.gov/fris/Home.aspx?ST=NC), where available; m-ain-stem reaches for which HEC-RAS crosssection data are unavailable will be assigned a representative cross section using best engineering judgment. Representative main-stem cross sections will be assigned based on the nearest available downstream main-stem cross section because a cross-section area generally increases from upstream to downstream. Tributary reaches for which cross-section data are unavailable will be assigned a representative tributary cross section based on the proximity to an available cross section and similar drainage area. After each reach is assigned the most appropriate cross section based on the location and drainage area, the discharge will be calculated for each reach by using length, slope, and cross-section data with the Manning’s equation shown in Equation 3-3. The channel slope (S ) for each reach will be calculated by dividing the difference between the maximum and minimum elevations by the reach length. =    2 2 3 11.486 Q A R Sn (3-3) where: 2 = Discharge (cfs) = Manning’s roughness coefficient = Cross-section area (squared feet [ft ]) = Hydraulic radius (ft) = Channel slope Q n A R S Manning’s roughness coefficients ()n have also been obtained from the North Carolina Flood Risk Information System. The values for the floodplain slope, channel slope, Manning’s roughness coefficient, and horizontal bank extension length will be set based on local topography and by using best engineering judgment, and the values can easily be adjusted during the calibration process. An F-table will be developed for each reach by calculating the surface area, volume, and discharge over a range of depths. The cross section can be extended 1,000 feet horizontally beyond each bank to allow the F-table to handle large storm flows. The floodplain slope will be assumed as 0.05. The volume and surface area will be calculated with the cross sections and stream segment lengths. The data used to calculate the elevation and slope for the model include the USGS 3D Elevation Program (https://www.usgs.gov/core-science-systems/ngp/3dep). 3.3 LAND SEGMENTATION Land-use, or land-cover, data is a critical factor in modeling watersheds because these data provide the detailed characterization of the potential pollutant sources entering the reaches as nonpoint-source contributions. Land-use distribution also has a major determining impact on the hydrologic response of the watershed. RSI-3348 22 2 This section describes how the proposed Pervious Land Segment (PERLND) and Impervious Land Segment (IMPLND) module-use categories were selected for explicit representation in the model application. The PERLND and IMPLND blocks of the UCI file contain most of the parameters that describe the way that water flows over and through the watershed. The objective of this task, therefore, will be to separate the watershed into unique land segments by using physical watershed characteristics to effectively represent the variability of hydrologic and water quality responses in the watershed. The primary watershed characteristics selected for the PERLND and IMPLND categorization include drainage patterns, meteorological variability, land cover, and soil properties. The municipal separate storm sewer system (MS4) areas will also be represented because of their link to permitting and water management. The NC DEQ has provided links to use for the MS4 areas in North Carolina. Final model land-cover characteristics will be selected based on the significance of their influence on the hydrologic processes and water quality constituents of interest, as well as the quality and availability of spatial data associated with the characteristics. 3.3.1 ELEVATION Topography provides elevation and slope values for the project area that are important to setting up HSPF because these values are needed for characterizing the landscape and land areas of the watershed. The flow accumulation and direction derived from elevation raster data are used to delineate subwatersheds. Average elevations and slopes are also calculated for each model subwatershed. The delineated subwatershed models are linked to the pervious or impervious lands that drain to the subwatersheds in the schematic block of the UCI file. Aggregating the subwatersheds into hydrozones based on meteorological variability will provide initial boundaries for the PERLNDs and IMPLNDs and allow for accurately representing the hydrologic processes while reducing computational demands. The procedures for determining the PERLND and IMPLND categories within each hydrozone are described in the following paragraphs. The 3D Elevation Program from the USGS has 10-m by 10-m elevation data across the United States available for download at (https://www.usgs.gov/core-science- systems/ngp/3dep). These 3D elevation data will be used to calculate the slope information for this model application. 3.3.2 LAND USE The NLCD 2019 land-cover layer (https://www.mrlc.gov) will be the primary layer used for the model land cover. Land covers will be aggregated/reclassified into a set of model land covers that will be used to develop the PERLND and IMPLND classifications within each subwatershed and hydrozone. Aggregated land-cover categories will be used to define the movement of water through the system (i.e., infiltration, surface runoff, and water losses from evaporation or transpiration) that is significantly affected by the land cover and its associated characteristics. 3.3.2.1 PERVIOUS AND IMPERVIOUS LAND CLASSIFICATION The number of operations (e.g., PERLND, IMPLND, RCHRES, PLTGEN, and COPY) allowed in one HSPF model application is limited; therefore, the categories represented in each state land-cover layer will be aggregated into relatively homogeneous model categories. Cropland will be separated using the National Agricultural Statistics Services Crop Data Layer (2021), as well as tile drainage status. Tile drainage status will be estimated as cropland with a slope of less than 4 percent calculated using USGS RSI-3348 23 2 3D elevation data (https://www.usgs.gov/core-science-systems/ngp/3dep) and hydrologic soil group (HSG) of A/D, B/D, C/D, or D from the Soil Survey Geographic Database (SSURGO, https://www.nrcs.usda.gov) Natural Resources Conservation Service (NRCS) Dataset [2022]. Figure 3-1 shows the general reclassification schemes for converting the land-cover classes to the model land- cover classes, and HSG categories are described in Table 3-1. HSG distributions by subwatershed can also be used as a basis for model parameterization related to infiltration and soil-moisture capacity values in the models, and the erodibility factor for each PERLND can be used to parameterize the erodibility factor of soils in the watershed. HSG percentages in each Neuse River Watershed project area are shown in Table 3-1. Figure 3-1. Land-Use Category Aggregation. Table 3-1. General Description of Hydrologic Soil Groups and Makeup of Watershed Areas Hydrologic Soil Group Abbreviated Description % Area A Sand; sandy loams with high-infiltration rates; well-drained soils with high transmission 20 B Silt loam or loam soils with moderate infiltration; moderately drained 14 C Sandy, clay loams; low-infiltration rates that impede water transmission 10 D Heavy soils, clay loams, silty, clay; low-infiltration rates that impede water transmission 8 AD A-Group soil, if drained 18 BD B-Group soil, if drained 20 CD C-Group soil, if drained 8 Unclassified No classification determined 3 RSI-3348 24 2 Feedlots located in the project area were identified using a layer downloaded from the NC DEQ Feature Server (https://services2.arcgis.com/kCu40SDxsCGcuUWO/ArcGIS/rest/services/Animal_Feed_ Operation_Permits_(View)/FeatureServer/0). Because this layer did not include poultry, poultry numbers were estimated by county with Census of Agriculture data. The feedlot data will be used to estimate fertilizer application to inform the calibration process throughout the watershed. Manure from feedlots will be assumed to have spread into the subwatersheds in which the feedlots are located within and adjacent to. The effective impervious area (EIA) is important to accurately represent in watershed models because of the EIA’s impact on the hydrologic processes that occur in urban environments. The term “effective” implies that the impervious region is directly connected to a local hydraulic conveyance system (e.g., gutter, curb drain, storm sewer, open channel, or river) and the resulting overland flow will not run onto pervious areas and, therefore, will not have the opportunity to infiltrate along the respective overland flow path before reaching a stream or waterbody. The average percent impervious on each developed model category will be calculated using the NLCD 2019 impervious layer. The data represent the total impervious area percentage, which will be used to determine the percent EIA by using Equation 3-4 from Sutherland [2000] to represent areas that are mostly storm sewered, with curb and gutter, and with residential rooftops connected to the MS4. ()=1.5EIA 0.1 TIA ,TIA 1 (3-4) Polluted stormwater runoff is commonly transported through MS4s before being discharged into local waterbodies. Certain MS4s are required to obtain National Pollutant Discharge Elimination System (NPDES) permits and develop stormwater management programs that describe stormwater control practices that will be implemented following permit requirements to minimize the discharge of pollutants from the storm sewer system [NPDES, 2023]. Representing regulated MS4s in the watershed in the HSPF model applications is important. GIS layers of the MS4 areas (i.e., polygons) were downloaded from the North Carolina Department of Transportation (NCDOT) (https://www.nconemap.gov/datasets/NCDOT::ncdot-city-boundaries/about) and NC DEQ (https://ncdenr.maps.arcgis.com/home/item.html?id=968f0ac608e44662bdb9a2d56350e605). MS4 areas will be represented in the model application schematic by using a separate mass link so that flow from those areas can be identified as separate from flow that originates in non-MS4 areas. 3.3.2.2 SEPTIC SYSTEMS A septic system falls under the category of on-site wastewater treatment systems (OWTS). OWTS are used by many households in the Neuse River Watershed. North Carolina has polygons that represent areas that are sewered (https://services.nconemap.gov/secure/rest/services/ NC1Map_Water_Sewer_2004/MapServer/2). Blockpop points that fall outside of the sewered areas provide the populations from the 2010 United States Census (https://www.census.gov/programs- surveys/decennial-census/data.html) that will be assumed to be on septic systems. OWTS are generally responsible for some pollutant loads to either the groundwater or tributaries. OWTS will be represented in the model application as a constant load and assumed to discharge at 50 gallons per person per day. The loading rates will initially be set at 40.4, 10.58, and 2.5 pounds per person per day for 5-day biochemical oxygen demand, total nitrogen, and total phosphorus, respectively [EPA,1980 and 1993]. The 5-day biochemical oxygen demand loads will be converted to carbonaceous biochemical oxygen RSI-3348 25 2 demand by using a factor of 1.2 for untreated waste [Thomann and Mueller, 1987]. Initial attenuation (i.e., pass-through) factors that represent septic-system efficiency will be set at 0.60, 0.77, and 0.14 for 5-day biochemical oxygen demand, total nitrogen, and total phosphorus, respectively [EPA,1980 and 1993; Vaudrey et al., 2016]. Soil attenuation will be represented as a function of simulated groundwater flow with less pass-through of pollutant loads occurring at lower flows, assuming more soil residence time results in greater pollutant degradation and transformation (e.g., denitrification). Initial failure rates and loading estimates will be source from the Chesapeake Assessment Scenario Tool (CAST; https://cast.chesapeakebay.net/About) from the southernmost Virginia area. RSI-3348 26 2 4.0 CALIBRATION AND VALIDATION 4.1 CALIBRATION AND VALIDATION TIME PERIODS Time-period selection for model calibration and validation depends on numerous factors, including the availability of data for model operations, land-use data for model setup, climate variability, and observed data for model-data comparisons. The principal time-series data that are needed for hydrologic and water quality calibration (i.e., meteorological, point-source, atmospheric deposition, observed flow, and water quality observations) indicate that a long-term simulation is possible at numerous streamflow gages throughout the Neuse River Watershed. Partial record periods, while not ideal, can still be used for consistency checks as part of the calibration and validation process. The date ranges for the calibration and validation periods include mixed wet and dry periods, as shown in Figure 4-1. Based on this consideration, the overall modeling period for hydrology and water quality is from 2003–2022, with 2003 being a “warm up” period and calibration beginning in 2004. Validation will likely be performed on the five most wet years and the five most dry years. Figure 4-1. Average Annual Precipitation for Modeled Watershed Areas from PRISM for Years 1990 to 2022. 4.2 HYDROLOGY CALIBRATION AND VALIDATION PROCEDURES AND COMPARISONS The Neuse River model application will be calibrated through an iterative process of making parameter changes, running the model, producing comparisons of simulated and observed values, and interpreting the results. This process will first occur for the hydrology portions of the model, followed by the water quality portions. The procedures have been well established during the past 35 years, as described in the application guide for HSPF [Donigian et al., 1984] and summarized by Donigian [2002]. The hydrology calibration process is greatly facilitated by using scripted processes in MATLAB. Calibrating HSPF to represent the hydrology of the Neuse River Watershed is an iterative trial-and-error process. Simulated results are compared with recorded data for the entire calibration period, including wet and dry conditions, to observe how well the simulation represents the hydrologic response under various climatic conditions. By iteratively adjusting specific calibration-parameter values within Driest Wettest Modeling Period RSI-3348 27 2 accepted and physically based ranges, the simulation results are changed until an acceptable comparison of simulation and recorded data is achieved [EPA, 2000]. The standard HSPF hydrologic calibration is divided into four phases: 1. Establish an annual water balance. This phase consists of comparing the total annual simulated and observed flows (in inches). The annual water balance is primarily governed by the input of rainfall and evaporation and the parameters for the lower zone nominal storage (LZSN), lower zone evapotranspiration parameter (LZETP), and infiltration index (INFILT). 2. Adjust low-flow/high-flow distribution. This step is generally performed by adjusting the groundwater or baseflow because the distribution between high and low flow is the easiest to identify in low-flow periods. Mean daily flow conditions are used, and the primary parameters involved are the INFILT, groundwater recession (AGWRC), and baseflow evapotranspiration index (BASETP). 3. Adjust storm-flow/hydrograph shape. The storm flow, which is compared in the form of short, timestep (1-hour) hydrographs, is largely composed of surface runoff and interflow. Adjustments are made with the upper zone storage (UZSN), interflow parameter (INTFW), interflow recession (IRC), and overland flow parameters (length of the overland flow plane [LSUR], Manning’s n [NSUR], and slope of the overland flow plane [SLSUR]). INFILT can also be used for minor adjustments. 4. Make seasonal adjustments. Differences in the simulated and observed total flow over each month and season are compared to see if runoff needs to be shifted from one month or season to another. These adjustments are generally accomplished by using seasonal (monthly variable) values for the parameters of vegetal interception (CEPSC), LZETP, and UZSN. Adjustments to variable groundwater recession (KVARY) and BASETP are also used. The procedures and parameter adjustments involved in these phases are more completely described in Donigian et al. [1984] and the HSPF hydrologic calibration expert system (HSPEXP) [Lumb et al., 1994; Duda et al., 2019]. The same model-data comparisons will be performed for the calibration and validation periods. The specific comparisons of simulated and observed values include: / Annual and monthly runoff volumes (inches) / Daily time series of flow (cubic feet per second [cfs]) / Storm-event periods (e.g., hourly values) (cfs) / Flow-frequency (flow-duration) curves (cfs) The water-balance components (input and simulated) are also reviewed. This effort involves displaying model results for individual land uses, as well as the entire watershed, for the following water-balance components: / Precipitation / Total runoff (sum of the following components): » Overland flow » Interflow » Baseflow / PEVT RSI-3348 28 2 / Total actual ET (sum of the following components): » Interception ET » Upper zone ET » Lower zone ET » Baseflow ET » Active-groundwater ET / Deep-groundwater recharge/losses Although observed values are not available for every water-balance component listed above, the average annual values must be consistent with expected values for the region, as impacted by the individual land-use categories. This consistency (or reality) check is separate, with data independent of the modeling (except for precipitation) to ensure that land-use categories and the overall water balance reflect the local conditions. The value ranges of the correlation coefficient (R) and R2 for assessing the model performance for daily and monthly flows are provided in Figure 4-2. The figure shows the range of values that may be appropriate for judging how well the model is performing based on the daily and monthly simulation results. As shown in Figure 4-2, the ranges for daily values are lower to reflect the difficulties in exactly duplicating the timing of flows given the uncertainties in the timing of model inputs, mainly precipitation. The general calibration and validation tolerances or targets that have been provided to model users as a part of HSPF training workshops over the past 20 years (e.g., Donigian [2000]) are listed in Table 4-1. The values in Table 4-1 attempt to provide general guidance in terms of the percent mean errors, or differences between simulated and observed values, so that users can determine what level of agreement or accuracy (i.e., very good, good, or fair) can be expected from the model applications. The target level of accuracy for this project will correspond in Table 4-1 to “Good” or “Very Good” results at more downstream main-stem calibration sites and “Fair” at more upstream tributary sites. Accuracy targets are highly dependent on the amount and quality of available data, and consequently, the targets will be finalized after the data gaps are analyzed. Figure 4-2. R and R2 Value Ranges for Model Performance. The caveats in the Table 4-1 notes indicate that the tolerance ranges should be applied to mean values and that individual events or observations may show larger differences and still be acceptable. The level of agreement to be expected also depends on numerous site- and application-specific conditions, including the data quality, purpose of the study, available resources, and available alternative assessment procedures that could meet the study objectives. RSI-3348 29 2 Table 4-1. General Calibration and Validation Targets or Tolerances for HSPF Applications [Donigian, 2000] Calibration Parameter Difference Between Simulated and Recorded Values (%) Very Good Good Fair Hydrology/Flow < 10 10–15 15–25 Stipulations: / Storm peaks may differ more than monthly and annual values / Quality detail of input and calibration data / Purpose of model application / Availability of alternative assessment procedures / Resource availability (i.e., time, money, and personnel) For any watershed modeling effort, the level of expected agreement is tempered by the complexities of the hydrologic system, quality of the available precipitation and flow data, and available information to characterize the watershed and quantify the human impacts on water-related activities. These tolerances are applied to comparisons of simulated and observed mean flows, annual runoff volumes, mean monthly and seasonal runoff volumes, and daily flow-duration curves. Larger deviations would be expected for individual storm events and flood peaks in both space and time. The values shown in Figure 4-2 were primarily derived from HSPF experience and past efforts on model performance criteria; however, the values do reflect common tolerances accepted by many modeling professionals. To provide a robust weight of statistical evidence, additional metrics, such as coefficient of model fit efficiency, root mean square error, seasonal differences, and low and high flow efficiency, will be reported. Additional evaluation criteria to be considered at primary gages are shown in Table 4-2. Table 4-2. Flow Calibration Criteria From Expert System for Calibration of HSPF [Lumb et al., 1994] Prediction Error Percent Difference Criteria (%) Error in Total Volume ±10 Error in Volume of 50% Lowest Flows ±10 Error in Volume of 10% Highest Flows ±15 Seasonal Volume Error (Summer) ±30 Seasonal Volume Error (Fall) ±30 Seasonal Volume Error (Winter) ±30 Seasonal Volume Error (Spring) ±30 Given the uncertain state of the art in model performance criteria, inherent errors in input and observed data, and approximate nature of model formulations, absolute criteria for watershed model acceptance or rejection are not generally considered appropriate by most modeling professionals. Most decision-makers, however, want a definitive answer to the question, “Is the model good enough for this evaluation?” Consequently, for the Neuse River modeling effort, the targets and tolerance ranges for RSI-3348 30 2 daily flows are proposed to correspond, at a minimum, to a fair to good agreement, and the ranges for monthly flows should correspond to good to very good agreement for calibration and validation at the primary calibration flow gages. Secondary calibration flow gages should ideally correspond to a fair agreement. Poor to fair ranges will be allowed for the tertiary sites because these sites are on smaller tributaries and usually have a much shorter representative dataset to work with. 4.3 WATER QUALITY CALIBRATION Water quality calibration is also completed through an iterative process of parameter adjustments and comparisons of simulated and observed values and is facilitated by using scripted processes in MATLAB. The model predictions are the integrated result of all of the assumptions used in developing the model input and representing the modeled sources and processes. Differences in model predictions and observations require the model user to reevaluate these assumptions for the estimated model input and parameters and consider the accuracy and uncertainty in the observations. Water quality monitoring sites will be mapped and assigned to each reach. Sites will be designated as primary, secondary, and tertiary (i.e., third priority). The most data-intensive, more downstream sites will be considered primary calibration sites. A calibration goal will be to keep the parameterization consistent throughout the project area to avoid curve fitting. Curve fitting is adjusting parameters reach by reach with the intent to force model results to follow the observed data curve without justification for why two neighboring reaches can exhibit such different behavior. Calibrating this way often causes many inconsistencies when using the model to define protection and restoration goals. In addition to completing the input development for the point-source, atmospheric deposition, and other contributions, the following steps will be performed at each of the calibration stations after the hydrologic calibration and validation: 1. Estimate all of the model parameters, which include land-use-specific accumulation and depletion/removal rates, wash-off rates, and subsurface concentrations. 2. Tabulate, analyze, and compare the simulated, annual, nonpoint loading rates with the expected range of nonpoint loadings from each land use (and each constituent) and adjust the loading parameters when necessary. 3. Calibrate instream water temperature, sediment, DO, and nutrients to the observed data. The primary calibration parameters involved in characterizing landscape-erosion processes are the coefficients and exponents from three equations that represent different soil detachment and removal processes [EPA, 2006]. Nonpoint sources of total ammonia and nitrate-nitrite will be simulated through accumulation and depletion/removal and a first-order wash-off rate from overland flow. Because of the affinity of orthophosphate to bind to sediments, orthophosphate will be simulated using a linear relationship with sediment washing off of the land. BOD will also be simulated using sediment-associated wash-off. Subsurface flow concentrations will be estimated on a monthly basis. Atmospheric depositions of nitrogen and ammonia will be applied to all of the land areas and contribute to the nonpoint-source load through the buildup/wash-off process. HSPF nonpoint loading rates (sometimes referred to as export coefficients) are highly variable, with values occasionally ranging to an order of magnitude depending on local and site conditions of soils, RSI-3348 31 2 slopes, topography, climate, and disturbance. As a quality check, the loadings from the Chesapeake Bay CAST system from the southernmost Virginia area will be compared to final loadings in the Neuse River HSPF model. The model simulates the instream and lake processes that contribute to sediment transport, algal growth, nutrient consumption, and DO dynamics. The sediment behavior for each size class will be investigated to ensure that the sediment dynamics reflect field observations. HSPF does not explicitly simulate streambank contribution dynamics; however, these processes will be implicitly included by allowing the streambed to contribute those loads. All of the required instream parameters will be specified for total ammonia, inorganic nitrogen, orthophosphate, and BOD. The processes in the instream portion of the models include BOD accumulation, storage, decay rates, benthic algal oxygen demand, settling rates, and reaeration rates. Atmospheric deposition onto water surfaces will be represented in the models as a direct input to the lakes and river systems. Biochemical reactions that affect DO will be represented in the model application. The overall sources considered for BOD and DO include point sources such as wastewater treatment facilities (WWTFs), nonpoint sources from the watershed, interflow, and active-groundwater flow. The instream calibration will begin with the temperature and sediment and then dissolved oxygen and nutrients. The dissolved oxygen and nutrient calibration will be conducted in tandem because these components depend on one another. The calibration requires developing time-series graphs to compare the simulated and observed water quality data. Instream water quality calibration will also include generating monthly boxplots, concentration-duration curves, and scatterplots of concentrations and corresponding flows. Hourly boxplots will be generated for temperature and dissolved oxygen to assess the diurnal variability. Sediment scour and deposition in the streambed for each reach over the period of simulation and the nutrient budget will also be evaluated. USGS groundwater concentration measurements will be used, when possible, for the water quality calibration. The parameters calibrated for sediment include buildup, wash-off, and instream scour and deposition. The parameters calibrated for nutrients include buildup, wash-off, instream cycling, algae growth, update, death, resuspension, and other parameters. Lake and impoundment water quality calibrations are often difficult in HSPF because the model simulates a completely mixed system (homogeneous waterbody with no vertical stratification). Large, deep headwater lakes with little inflow/outflow can produce an overestimation/accumulation of nitrogen or phosphorus and low chlorophyll a concentrations. The inherent variability in depth, surface area, volume, and residence time between lakes in a model application causes each lake to behave differently, which makes developing a standard approach or solution difficult. To address these issues and achieve a dynamic, steady-state system, instream parameters for lakes and ponds are generally very different compared to reaches [AQUA TERRA Consultants, 2015]. Bias as a result of lake stratification is reduced by calibrating to the observed values taken from the surface to 1 m in depth. The main goal of watershed water quality calibration is to obtain acceptable agreement of observed and simulated concentrations (i.e., within defined criteria or targets) while keeping the instream water quality parameters within physically realistic bounds and the nonpoint loading rates within the expected ranges RSI-3348 32 2 from the literature. The general water quality calibration targets or tolerances for HSPF applications are shown in Table 4-2. The calibration should be accomplished while maintaining consistent parameters in each land-use category throughout the watershed, when possible. Table 4-3. General Calibration Targets or Tolerances for HSPF Applications [Donigian, 2000] Calibration Parameter Difference Between Simulated and Recorded Values (%) Very Good Good Fair Sediment < 20 20–30 30–45 Water Temperature < 7 8–12 13–18 Water Quality/Nutrients < 15 15–25 25–35 Pesticides/Toxics < 20 20–30 30–40 Stipulations: / Storm peaks may differ more than monthly and annual values / Quality detail of input and calibration data / Purpose of model application / Availability of alternative assessment procedures / Resource availability (i.e., time, money, personnel) 4.4 SENSITIVITY ANALYSIS Traditional sensitivity analysis involves comparing the relative differences in the results by dividing the percent change in a variable response by the percent change in a calibration parameter, which provides an understanding for how parameters impact the physical and biological processes represented in the HSPF model. Although this analysis is valuable and inherent during the calibration process, several thousands of parameter-variable options exist for evaluation, and interpreting the results requires a thorough understanding of HSPF mechanics and terminology. A more useful approach to understanding model sensitivity involves a detailed evaluation of source contributions using real-world terminology. By using the HSPF model application results for each day, source, basin, and constituent of concern, source contributions at all of the locations for selected time periods and/or flow conditions can be calculated. The source contribution analysis will be included for the outlet of each model application (i.e., all of the sources contributing to a main-stem river endpoint) as part of the final model deliverables. A more detailed source analysis can be done with the SAM tool, which will be included as part of the HSPF and SAM trainings. RSI-3348 33 2 5.0 DATA MANAGEMENT Two data types will be used to support the Neuse River HSPF modeling project: GIS and time-series data. The data types must change format as they are integrated into an HSPF model and are thus subject to possible errors. RESPEC will adhere to electronic data acquisition protocols that were developed for previous HSPF applications. These protocols ensure that quality assurance considerations related to preventing, detecting, and correcting electronic data manipulation errors are properly addressed. RESPEC will also adhere to the protocols for data acceptance criteria described in the Neuse River Watershed Model Quality Assurance Project Plan (QAPP) [Duda and Kirby, 2023]. RESPEC will maintain a copy of the project files on the network for a minimum of 5 years following the completion of the project. Consistent data management procedures will be used during the preprocessing, model calibration, and postprocessing stages of the project. All of the data and information collected and generated during this project will be stored in a project folder on RESPEC’s network. Data processing will be completed using a combination of ArcGIS, MATLAB, Python, and the SARA Time-Series Utility. RESPEC modelers will be responsible for adhering to and documenting data management practices that ensure the quality of the data that are downloaded and/or manipulated. Original data sources will be documented to identify the website or contact person that provided the data, data query parameters, and data request correspondence. Original, unaltered copies of all data sources used in the project will be retained in the project folder on RESPEC’s network. Metadata will be included with spatial datasets. The SARA Time- Series Utility will be used to access the WDM files that will be used to store model-input data such as meteorological, point-source, atmospheric deposition, and other time-series data. GIS data will be used in a geodatabase feature-class format. The projection of all of the GIS data will be consistent. When new GIS data are added to a feature class, ArcPro automatically projects the data to match the projection of the feature class. Example metadata standards for NC OneMap are available online (https://www.nconemap.gov/pages/metadata). Model inputs, including meteorological data, point-source data, and surface-water withdrawals, will be stored in a WDM file during the calibration process. Model outputs at calibration gages will be stored in a set of binary (HBN) files. The HBN and WDM files can be accessed using the SARA Time-Series Utility, which can be downloaded online (https://www.respec.com/product/modeling-optimization/sara- timeseries-utility). SAM files, which are downloaded to the user’s local drive, will also be used as storage, and the SAM program will have the capability to allow the user to extract the model data. RSI-3348 34 2 6.0 REFERENCES AQUA TERRA Consultants, 2015. NPS Target Loading Rates for Minnesota, prepared by AQUA TERRA Consultants, Mountain View, CA, for the Minnesota Pollution Control Agency, St. Paul, MN. Bicknell, B. R.; J. C. Imhoff; J. L. Kittle, Jr.; T. H. Jobes; and A. S. Donigian, Jr., 2005. Hydrological Simulation Program–Fortran (HSPF), User's Manual for Release 12.2, prepared by the U.S. Environmental Protection Agency, Athens, GA. Donigian, Jr., A. S., J. C. Imhoff, B. R. Bicknell, and J. L. Kittle, Jr., 1984. Application Guide for Hydrological Simulation Program - Fortran (HSPF), EPA-600/3-84-065, prepared by Environmental Research Laboratory, Office of Research and Development, U.S. Environmental Protection Agency, Athens, GA. Donigian, Jr., A. S., B. R. Bicknell, J. C. Imhoff, J. L. Kittle, Jr., T. H. Jobes, and P. B. Duda, 2018. HSPF Version 12.5 User’s Manual, prepared by RESPEC, Inc., Mountain View, CA. Duda, P. B., A. Mishra, and A. S. Donigian, Jr., 2019.HSPEXP+ User’s Manual, Version 3.0, RSI-2770, prepared by RESPEC, Rapid City, SD, for the Minnesota Pollution Control Agency, St. Paul, MN. Donigian, Jr., A. S., 2000. HSPF Training Workshop Handbook and CD, Lecture #19: Calibration and Verification Issues, Slide #L19-22, prepared by the U.S. Environmental Protection Agency Headquarters, Washington Information Center, for the U.S. Environmental Protection Agency, Office of Water, Office of Science and Technology, Washington, DC. Donigian, Jr., A. S., 2002. “Watershed Model Calibration and Validation: The HSPF Experience,” Proceedings, Water Environment Federation National Total Maximum Daily Load Science and Policy Conference, Phoenix, AZ, November 13−16. Chapra, S. C., 1997. Surface Water-Quality Modeling, Waveland Press, Inc., Long Grove, IL. Duda, P. B. and C. M. Kirby, 2023. Neuse River Watershed Model Quality Assurance Project Plan, QAPP- 21, prepared by RESPEC, Rapid City, SD, for the North Carolina Department of Environmental Quality, Raleigh, NC. Gardner, C. B., P. J. Zarriello, G. E. Granato, J. P. Masterson, D. A. Walter, A. M. Waite, and P. E. Church, 2011. Simulated Effects of Water Withdrawals and Land-Use Changes on Streamflows and Groundwater Levels in the Pawcatuck River Basin, Southwestern Rhode Island and Southeastern Connecticut, Scientific Investigations Report 2009–5127, prepared by the U.S. Department of the Interior, U.S. Geological Survey, Reston, VA; U.S. Department of Agriculture, Natural Resources Conservation Service, Washington, DC; and Rhode Island Water Resources Board, Providence, RI. Hu, H. L., H. M. Chen, N. P. Nikolaidis, D. R. Miller, and X. Yang, 1998. “Estimation of Nutrient Atmospheric Deposition to Long Island Sound,” Water, Air, and Soil Pollution, Vol. 105, pp. 521–538. Kenner, S. J., 2023. Selection of HSPF Watershed Model for the Neuse River. RSI(RCO)- W0392.22001/3-23/2, memorandum from S. Kenner, RESPEC, Rapid City, SD, to P. Behm, North Carolina Department of Environmental Quality, Raleigh, NC, March 1. Koelliker, Y., L. A. Totten, C. L. Gigliotti, J. H. Offenberg, J. R. Reinfelder, Y. Zhuang, and S. J. Eisenreich, 2004. “Atmospheric Wet Deposition of Total Phosphorus in New Jersey,” Water, Air, and Soil Pollution, Vol. 154, pp. 139–150. RSI-3348 35 2 Kohler, M. A., T. J. Nordenson, and W. E. Fox, 1955. Evaporation From Pans and Lakes, Research Paper 38, prepared by the U.S. Department of Commerce, Weather Bureau, Washington DC. Lumb, A. M.; R. B. McCammon; and J. L. Kittle, Jr., 1994. Users Manual for an Expert System (HSPEXP) for Calibration of the Hydrological Simulation Program-FORTRAN, U.S. Geological Survey Water Resources Investigations Report 94-4168, prepared by the U.S. Geological Survey, Reston, VA. Menne, M. J., I. Durre, B. Korzeniewski, S. McNeal, K. Thomas, X. Yin, S. Anthony, R. Ray, R. S. Vose, B. E. Gleason, and T. G. Houston, 2022. “Global Historical Climatology Network - Daily (GHCN-Daily), Version 3.22,” noaa.gov, accessed November 20, 2022, from ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/daily National Pollutant Discharge Elimination System, 2023. “Stormwater Discharges From Municipal Sources,” epa.gov, accessed April 24, 2023, from https://www.epa.gov/npdes/stormwater-discharges- municipal-sources Natural Resources Conservation Service, 2022. “Web Soil Survey (WSS),” usda.gov, accessed November 30, 2022, from https://websoilsurvey.nrcs.usda.gov/ Penman, H. L., 1948. “Natural Evaporation From Open Water, Bare Soil, and Grass,” Royal Society, Vol. 193, No. 1032, pp. 120–145. Stull, R., 2017. Practical Meteorology: An Algebra-Based Survey of Atmospheric Science, Version 1.02b, published by the University of British Columbia, Vancouver, BC, Canada. Sutherland, R. C., 2000. “Methods for Estimating the Effective Impervious Area of the Urban Watersheds,” Technical Note #28, Watershed Protection Techniques, Vol. 2, No. 1, pp. 282–284. Thompson, E. S., 1976. “Computation of Solar Radiation From Sky Cover,” Water Resources Research, Vol. 12, No. 5, pp. 859–865. Thomann, R. V. and J. A. Mueller, 1987. Principles of Surface Water Quality Modeling and Control, Harper-Collins, New York, NY. U.S. Bureau of Reclamation, 1987. Design of Small Dams, 3rd Edition, published by the U.S. Government Office of Printing, Washington, DC. U.S. Environmental Protection Agency, 1980. Design Manual, Onsite Wastewater Treatment and Disposal Systems, EPA-625/1-80-012, prepared by the U.S. Environmental Protection Agency, Office of Water Program Operations, Washington, DC. U.S. Environmental Protection Agency, 1993. Guidance Specifying Management Measures for Sources of Nonpoint Pollution in Coastal Waters, EPA-840-B-92-002, prepared by the U.S. Environmental Protection Agency, Office of Water Program Operations, Washington, DC. U.S. Environmental Protection Agency, 2000. EPA BASINS Technical Note 6: Estimating Hydrology and Hydraulic Parameters for HSPF, EPA-823-R00-012, prepared by U.S. Environmental Protection Agency, Office of Water Program Operations, Washington, DC. U.S. Environmental Protection Agency, 2006. EPA BASINS Technical Note 8: Sediment Parameter and Calibration Guidance for HSPF, prepared by U.S. Environmental Protection Agency, Office of Water Program Operations, Washington, DC. RSI-3348 36 2 Vaudrey, J. M. P, C. Yarish, J. K. Kim, C. Pickerell, L. Brousseau, J. Eddings, M. Sautkulis, 2016. Connecticut Sea Grant Project Report: Comparative Analysis and Model Development for Determining the Susceptibility to Eutrophication of Long Island Sound Embayments, Project number R/CE-34-CTNY, prepared for the Connecticut Sea Grant, Groton, CT; New York Sea Grant, Stony Brook, NY; and the Long Island Sound Study, Stamford, CT. World Meteorological Organization, 2014. Guide to Meteorological Instruments and Methods of Observation, WMO-No. 8, prepared by the World Meteorological Organization, Geneva, Switzerland. Yang, X., D. R. Miller, X. Xu, L. H. Yang, H. M. Chen, and N. P. Nikolaidis, 1996. “Spatial and Temporal Variations of Atmospheric Deposition in Interior and Coastal Connecticut,” Atmospheric Environment, Vol. 30, No. 22, pp. 3801–3810.