USVI Public Records

Document text

APPENDIX D: St. Thomas East End Reserve, St. Thomas, USVI Regional Model Development Report JUNE 2022 PREPARED FOR: US Environmental Protection Agency Region 2 290 Broadway New York, NY 10007-1866 PREPARED BY: RTI International 3040 Cornwallis Road Research Triangle Park, NC 27709 Paradigm Environmental 4021 University Drive Suite 203 Fairfax, VA 22030 Appendix D – Modeling Report ii June 2022 THIS PAGE INTENTIONALLY LEFT BLANK Appendix D - Modeling Report June 2022 iii CONTENTS Acronyms and Abbreviations ....................................................................................................... vi 1 Problem Definition and Management objectives ..................................................................... 1 1.1 Background ................................................................................................................... 1 2 Conceptual Model: key processes & variables ......................................................................... 2 2.1 The Model Development Cycle ...................................................................................... 2 3 Technical approach ............................................................................................................... 4 3.1 Watershed Modeling System .......................................................................................... 4 3.2 TMDL Mixing Model .................................................................................................... 5 3.2.1 Conceptual Model .................................................................................................. 6 4 Model Development ............................................................................................................. 8 4.1 Hydrologic Response Units ............................................................................................ 8 4.1.1 Land Use/Land Cover ............................................................................................ 8 4.1.2 Soils Data ............................................................................................................. 10 4.1.3 Slope .................................................................................................................... 12 4.1.4 Other HRU Refinements ....................................................................................... 14 4.2 Source Characteristics .................................................................................................. 19 4.3 Data Availability and Quality ....................................................................................... 22 4.3.1 Meteorological Data ............................................................................................. 23 4.3.2 Potential Evapotranspiration ................................................................................. 27 4.4 Data Gaps ................................................................................................................... 28 4.4.1 Hydrological Data Gaps ........................................................................................ 28 4.4.2 Water Quality Data Gaps ...................................................................................... 28 4.5 Time Frame of Simulation............................................................................................ 30 5 Model Calibration ............................................................................................................... 30 5.1 Water Balance Calibration ........................................................................................... 31 5.2 Model Parameters ........................................................................................................ 33 5.3 Water Quality Calibration ............................................................................................ 34 6 Model evaluation and acceptance ........................................................................................ 36 6.1 Calibration Example: Jolly Hill Gut, St. Croix .............................................................. 36 6.2 Calibration Example: Turpentine Run, St. Thomas ....................................................... 39 6.3 Validation Example: Turpentine Run, St. Thomas ......................................................... 42 6.4 Water Quality Model ................................................................................................... 43 7 References .......................................................................................................................... 47 Appendix D – Modeling Report iv June 2022 Figures 2-1. Conceptual schematic of a model development cycle. ............................................................ 3 3-1. LSPC water budget schematic. ............................................................................................. 5 3-2. Global tidal variation (Source: Dr. Richard Ray/Space Geodesy branch, NASA/GSFC, NASA: http://svs.gsfc.nasa.gov/stories/topex/images/TidalPatterns_hires.tif) ................................ 7 3-3. Conceptual model for estimating critical-condition mixing volume for each water body. ......... 7 4-1. NOAA 2015 C-CAP land use/land cover (top panel) and NOAA pervious/impervious cover with TIGER Census roads layer (bottom panel) for Saint Thomas, USVI. ............................. 9 4-2. STATSGO/SSURGO hydrologic soil groups for St. Thomas, USVI.................................... 11 4-3. National Elevation Dataset coverage for St. Thomas, USVI. ............................................... 12 4-4. National Elevation Dataset slope coverage for St. Thomas, USVI. ....................................... 13 4-5. Standard United States Department of Agriculture (USDA) Soil Triangle with Hydrologic Soil Group mapping. ............................................................................................................... 17 4-6. STATSGO/SSURGO K-Factor for St. Thomas, USVI. ...................................................... 17 4-7. NCDC meteorological gaging stations on St. Thomas, USVI. .............................................. 23 4-8. Total Precipitation at CHARLOTTE AMALIE CYRIL (671640), 1991-2015. ..................... 26 4-9. Total Precipitation at CHRISTIANSTED AP (671624), 1991-2015. .................................... 26 4-10. Regional NCDC rainfall totals vs. gage elevation. ............................................................... 27 5-1. Process for model calibration to minimize propagation of uncertainty. ................................. 30 5-2. Model parameterization and calibration sequence for land hydrology. ................................. 32 6-1. Temporal data quality summary for model calibration at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). ......................................................................... 37 6-2. Daily modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). ............................................................................... 37 6-3. Monthly modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). ............................................................................... 38 6-4. Seasonal average modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). ......................................................................... 38 6-5. Seasonal interquartile modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). .............................................................. 39 6-6. Temporal data quality summary for model calibration at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000)...................................................................... 40 6-7. Daily modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). ........................................................................... 41 6-8. Monthly modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). ........................................................................... 41 6-9. Seasonal average modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000)...................................................................... 42 6-10. Seasonal interquartile modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). ......................................................... 42 6-11. Modeled unit-area sediment (lbs/acre) for all watersheds..................................................... 44 6-12. Modeled unit-area BOD5 (lbs/acre) for all watersheds. ........................................................ 44 6-13. Modeled unit-area ENTERO (Billion/acre) for all watersheds. ............................................ 45 6-14. Modeled unit-area TN (lbs/acre) for all watersheds. ............................................................ 45 6-15. Modeled unit-area TP (lbs/acre) for all watersheds. ............................................................. 46 Appendix D – Modeling Report June 2022 v Tables 1-1. List of waterbodies included in this project and corresponding 303(d) listed constituents ......... 1 4-1. Composite land use categories and model groups (Source: 2015 NOAA C-CAP) ................. 10 4-2. Land use and soils overlay (Source Data: 2015 NOAA C-CAP and STATSGO/SURGO Soils) ............................................................................................................................ 11 4-3. Land use and slope overlay (Source Data: 2015 NOAA C-CAP and NED DEM slopes) ...... 13 4-4. Land use change matrix (Source Data: 2007 and 2015 NOAA C-CAP land use) .................. 15 4-5. K-Factor vs. hydrologic soil group and associated soil textures (Source Data: NRCS 2016a, 2016b) .............................................................................................................................. 16 4-6. Mapped Impervious Area (MIA) vs. Effective Impervious Area (EIA) summarized by TMDL watershed ......................................................................................................................... 18 4-7. Summary of regional and national unit-area loads by land use (LU) from various literature sources ............................................................................................................................ 20 4-8. Data quantity and quality inventory for model development ................................................ 22 4-9. Summary of Virgin Islands NCDC rainfall data compiled and processed for modeling ......... 25 4-10. Estimated stratification of modeled ET by land cover type ................................................... 28 4-11. Quantity and quality of USGS streamflow records by water year (Source Data: USGS) ........ 29 5-1. Upstream drainage area HRU characteristic for each USGS gage ........................................ 31 5-2. Typical performance targets for LSPC hydrology simulation (modeled vs. observed) ............ 33 5-3. Summary of calibrated hydrology parameters ..................................................................... 33 5-4. Legend for geometric properties table ................................................................................. 34 5-5. Summary of measured geometric properties and calibrated critical mixing ratio and volume for each Assessment Unit ....................................................................................................... 35 6-1. Model performance at Jolly Hill Gut for good-and-estimated rainfall vs. flow data quality intervals ............................................................................................................................ 36 6-2. Model performance for Turpentine Run at Mt. Zion, St. Thomas, for good-and-estimated rainfall vs. flow data quality intervals ................................................................................. 40 6-3. Comparison of the mean annual and 95th percentile flow between StreamStats and model for Turpentine Run at Mt. Zion, St. Thomas ........................................................................... 43 6-4. Critical condition month for each TMDL waterbody and pollutant under existing conditions43 Appendix D – Modeling Report vi June 2022 ACRONYMS AND ABBREVIATIONS BOD Biochemical Oxygen Demand BOD5 Biochemical Oxygen Demand (5-day) C-CAP Coastal Change Analysis Program DCIA Directly Connected Impervious Area DEM Digital Elevation Model DO Dissolved Oxygen DPNR Department of Planning and Natural Resources (U.S. Virgin Islands) ENTERO Enterococcus ET Evapotranspiration GIS Geographic Information System HRU Hydrologic Response Unit HSPF Hydrologic Simulation Program--FORTRAN ICIS Integrated Compliance Information System LiDAR Light Detecting and Ranging LSPC Loading Simulation Program C++ MIA Mapped Impervious Area NASA National Aeronautics and Space Administration NCDC National Climatic Data Center NED National Elevation Dataset NHDPlus National Hydrography Dataset Plus NOAA National Oceanic and Atmospheric Administration NPDES National Pollutant Discharge Elimination System QAPP Quality Assurance Project Plan QCLCD Quality-Controlled Local Climatological Data QGIS Quantum Geographic Information System (software) SOD Sediment Oxygen Demand SSURGO Soil Survey Geographic database STEER St. Thomas East End Reserve STORET EPAs Storage and Retrieval System STATSGO State Geographic Database STXEEMP St. Croix East End Marine Park SUSTAIN System for Urban Stormwater Treatment and Analysis Integration Model TMDL Total Maximum Daily Load TN Total Nitrogen TP Total Phosphorus TSS Total Suspended Sediment USDA United States Department of Agriculture USEPA U.S. Environmental Protection Agency USFS United States Forest Service USGS United States Geological Survey USVI United States Virgin Islands WAPA Water and Planning Authority WTM Watershed Treatment Model WQS Water Quality Standards Appendix D - Modeling Report June 2022 1 1 PROBLEM DEFINITION AND MANAGEMENT OBJECTIVES This document provides a summary of the methods and modeling approaches utilized in support of the Total Maximum Daily Load (TMDL) calculations for the waterbodies listed in Table 1-1. Details on the data used in developing the modeling extent and parameterization is presented and model calibration and verification results are included. This document is intended to be used as a supplement to the TMDL reports and the reader is encouraged to read each individual TMDL report for descriptions of the watersheds, project objectives, and results. This document explains the major steps, decisions, and assumptions made in the watershed modeling process. The STEER TMDLs were initially developed in parallel with other TMDLs on St. John and St. Croix. A regional modeling approach was used in which a model was set up and calibrated to provide consistent coverage of all three islands. While this modeling report focuses on STEER watersheds and waterbodies specifically, there are some elements of the model development process that are presented for watersheds outside of STEER (e.g., climate data in Section 4.3.1, calibration discussion in Section 5) due to the availability, or lack of, STEER specific data. 1.1 Background Table 1-1 summarizes the waterbodies included on the United States Virgin Islands (USVI) 303(d) list of impaired waters requiring TMDLs. A variety of pollutants and impairment types are included, some of which can be directly related to water quality concentrations and some of which must be related via a surrogate target (i.e., Biochemical Oxygen Demand (BOD5), Total Nitrogen (TN) and Total Phosphorous (TP) loadings impacting Dissolved Oxygen (DO); Total Suspended Sediment (TSS) impacting Turbidity; and Enterococcus (ENTERO) as a surrogate for general indicator bacteria). The watershed model developed to support these TMDLs was set up based on our understanding of these systems, the major sources, and how pollutants are delivered to the waterbodies. Specific information on each of these impairments can be found within the TMDL reports. Table 1-1. List of waterbodies included in this project and corresponding 303(d) listed constituents Watershed Waterbody ID Waterbody Name Dissolved Oxygen Enterococcus Bacteria Fecal Coliform pH Turbidity St. Thomas East End Reserve (STEER) VI-STT-25 Great Bay  VI-STT-28 Cowpet Bay  VI-STT-31 Nazareth Bay  VI-STT-33 Benner Bay VI-STT-34 Benner Bay Lagoon Marina   VI-STT-35 Mangrove Lagoon   Appendix D – Modeling Report 2 June 2022 2 CONCEPTUAL MODEL: KEY PROCESSES & VARIABLES The primary objective of this study is to support TMDL calculations for implementation planning. Using all available data, this modeling approach quantifies flow and loads from anthropogenic and natural sources. The watershed model provides a high-resolution assessment of sources and loading potential. For TMDL calculation, the linkage analysis represents the essence of critical conditions associated with the TMDL impairment within the waterbody. This modeling approach uses a top- down weight of evidence approach to correlate the best-available information and corroborate findings across different layers of the analysis. 2.1 The Model Development Cycle For TMDL analysis, the model development process can be a good platform for gaining valuable information and insight about a natural system. If well-designed, the model development process is an iterative and adaptive cycle that improves understanding of the natural system over time as better information becomes available. Ultimately, a model can inform future data acquisition efforts and management decisions by highlighting factors that have the most impact on the behavior of a natural system. Figure 2-1 is a conceptual schematic of a model development cycle. A well-designed model development cycle is conceptually circular as opposed to linear. That cycle can be summarized in six interrelated steps: 1. Assess Available Data: these data are used for source characterization, trends analysis, and defining modeling objectives 2. Delineate Project Extent: which refers to model segmentation and discretization 3. Set Boundary Conditions: including quality-controlled spatial and temporal model inputs 4. Represent Processes: refers to calibration of model rates and constants to mimic observed physical processes of the natural system 5. Confirm Responses: refers to validation of model processes over space and time to assess if the model is a robust predictive tool. 6. Assess Data Gaps: Sometimes the rigidity of modeled responses can highlight unrepresented physical processes in the natural system. Those data gaps sometimes provide a sound basis for further data collection efforts to refine the model, which cycles back to Step 1. Appendix D – Modeling Report June 2022 3 Figure 2-1. Conceptual schematic of a model development cycle. Appendix D – Modeling Report 4 June 2022 3 TECHNICAL APPROACH This section describes the models used as the scientific basis to establish current conditions and evaluate water quality-based controls to reduce pollution from both point and nonpoint sources and to restore and maintain the quality of the state's water resources (USEPA, 1991). 3.1 Watershed Modeling System The hydrologic and water quality model applied for this TMDL was the Loading Simulation Program in C++ (LSPC), a watershed modeling system that includes Hydrologic Simulation Program– FORTRAN (HSPF) algorithms for simulating watershed hydrology, erosion, water quality processes, and in-stream fate and transport processes. LSPC integrates Geographic Information System (GIS) outputs, comprehensive data storage and management capabilities, the original HSPF algorithms, and a data analysis/post-processing system into a convenient PC-based Windows environment. The algorithms of LSPC are identical to a subset of those in the HSPF model with selected additions, such as algorithms to address land use change over time. LSPC is a public domain watershed model originally made available through United States Environmental Protection Agency’s (USEPA’s) Office of Research and Development in Athens, Georgia as a component of USEPA’s National TMDL Toolbox. Some of the most recent advancements and applications of LSPC are in the Los Angeles, California Region (http://dpw.lacounty.gov/wmd/wmms/). Simulation of upland loading and transport of sediment and nutrients was conducted to develop estimates of the relative contribution of pollutants by source to support developing land-based mitigation strategies. Three previous modeling efforts were reviewed and referenced for this study: 1. DO TMDL for Salt River Bay, St. Croix (DPNR, 2004). Continuous-simulation modeling using the LSPC watershed model. 2. Watershed Characterization and Planning for Pathogen Source Reduction in the USVI (Cadmus, 2011). Continuous-simulation modeling using the System for Urban Stormwater Treatment and Analysis INtegration (SUSTAIN). 3. St. Croix East End Watersheds Management Plan (Horsley Witten Group, 2011). Planning- level spreadsheet modeling using the Watershed Treatment Model (WTM). The first two studies are process-based continuous simulation approaches that are capable of providing predicted flow and loading for a wide range of conditions that vary over space and time. Findings from those two studies provided directly applicable reference material for parameterizing and calibrating parameters associated with hydrologic and water quality processes, while the third provided locally derived model coefficients and long-term estimated source loads for benchmark comparison of simulated model results. The hydrologic (water budget) model in LSPC simulates the stages of the water cycle for each Hydrologic Response Unit (HRU) segment, from rainfall on constructed landscapes, vegetation, and soil. Variations in interception storage, surface conditions, soil type, slope, length of overland flow, and other physical parameters influence the fate of water as it moves through the various layers of land, as applicable. Different HRUs will respond in different ways depending on how each is parameterized. Within the LSPC Model environment water exits the land segment and is routed to stream segments by up to three different flow pathways: surface runoff, interflow outflow, and/or groundwater outflow. Figure 3-1 presents a schematic of model inputs; storage compartments; transfers, losses, and interactions; and parameters (rates and constants). Appendix D – Modeling Report June 2022 5 Figure 3-1. LSPC water budget schematic. 3.2 TMDL Mixing Model The listed TMDL segments are coastal estuarine waterbodies. Observed data for water quality were collected at different stations within the waterbodies, but not upstream in the watershed. The spatial and temporal range of observed data represent different mixing conditions in each of the waterbodies, with some waterbodies providing a higher dilution potential than others. The objective of the TMDL is to capture the critical condition associated with the impairment and express representative daily load and waste load allocations for implementation planning. A review of locally developed TMDLs showed that several of the DO TMDLs include Sediment Oxygen Demand (SOD) loading as a significant source of the impairment (St. Thomas Harbor, Vessup Bay and Red Hook Bay, Mangrove Lagoon, Benner Bay Lagoon Marina, and Benner Bay, St. Thomas; North Shore, and Salt River Bay, St. Croix). SOD represents an accumulation of loads from contributing sources at one point in time that affects water quality conditions at a later point in time. USVI is in a region of the planet where tidally influenced mixing is low. Upwelling and wind-driven storm mixing are the primary drivers of water exchange within the estuarine waterbodies. Building upon insights gained from previous regional TMDL modeling efforts, along with trends analysis to characterize critical conditions associated with water quality impairments, a simple mixing model was applied to represent critical conditions and calculate daily loads from contributing sources for the TMDL. Appendix D – Modeling Report 6 June 2022 3.2.1 Conceptual Model At the core of all hydrodynamic and water quality models are the fundamental concepts of continuity and mass balance. Previous modeling efforts have demonstrated the complex interactions associated with dilution and mixing, reactions and transformations. For mass balance of conservative substances in estuarine conditions, the continuity equation to characterize mixing takes the following form: 𝑄𝑖𝑛𝑆0 + 𝑄𝑓𝑆𝑓= 𝑄𝑜𝑢𝑡𝑆1 where: Qin = Oceanic Inflow S0 = Oceanic salinity Qf = Freshwater inflow Sf = Freshwater salinity (generally, Sf = 0) Qout = Outflow to the ocean S1 = Outflow salinity Commonly used indicators of estuary health include residence time and exposure time. Residence time represents the time it takes for a water drop to leave the estuary, while exposure time represents the time a drop spends in the estuary until it never returns (which accounts for water that re-enters the estuary during ebb/flow cycles). A well-flushed estuary generally has better water quality than a poorly flushed estuary. Residence time (T) is expressed as: 𝑇= 𝑉𝑜𝑙 𝑄𝑓 (1 − 𝑆1 𝑆0 − 𝑆𝑓 𝑆0 ) where: Vol = Volume of estuary within estuary Qf = Freshwater inflow S0 = Oceanic salinity Sf = Freshwater salinity (generally, Sf = 0) S1 = Outflow salinity Exposure time (t) is expressed as: 𝑡= 𝑉𝑜𝑙̅̅̅̅̅ × 𝑇𝑡𝑖𝑑𝑒 (1 −𝑟)𝑉𝑇𝑃 where: 𝑉𝑜𝑙̅̅̅̅̅ = Mean volume of estuary Ttide = Tidal period r = Return coefficient (r=1, all re-enters; r=0, none re-enters) VTP = Mean tidal prism volume As shown in Figure 3-2, the Caribbean is located at an amphidromic point on the Earth, meaning that it experiences high tide all the time; so effectively, there is negligible tide because the water level is nearly constant. Amphidromic points are the dark blue areas shown in Figure 3-2 where the white lines converge. Because the tidal period and volume of tide are less applicable to the Caribbean, the major factors affecting exposure time are the volume of water in estuary (Figure 3-3) and the return coefficient (which is a value between 0 and 1). Some physical factors that influence both residence time and exposure time are: 1. The width (or narrowness) of the estuary opening relative to its shoreline perimeter 2. The size of the watershed/drainage area (i.e., freshwater inflow volume) 3. The orientation of the estuary’s opening relative to upwelling or wind-driven ocean currents Appendix D – Modeling Report June 2022 7 4. Bathymetric features or obstructions (e.g. sills, corals, other islands) Figure 3-2. Global tidal variation (Source: Dr. Richard Ray/Space Geodesy branch, NASA/GSFC, NASA: http://svs.gsfc.nasa.gov/stories/topex/images/TidalPatterns_hires.tif) After the existing condition was established, load allocations were applied to individual sources until modeled concentrations achieved target concentrations for the TMDL. TMDL calculation details, as well as summary tables with TMDL information for each listed waterbody and constituent, are included in the respective TMDL documents for each of the watersheds. Figure 3-3. Conceptual model for estimating critical-condition mixing volume for each water body. Appendix D – Modeling Report 8 June 2022 4 MODEL DEVELOPMENT 4.1 Hydrologic Response Units An HRU is the smallest spatial unit in the model. Each HRU represents areas of similar physical characteristics attributable to certain processes. Spatial or geological characteristics such as land use, soils, and slopes are typically used to define HRUs. This section describes each of the component layers and the overlay process used to derive the HRUs. 4.1.1 Land Use/Land Cover Land use/land cover data are the primary layers that form the basis for HRUs. The primary sources of land use/land cover for this effort were the National Oceanic and Atmospheric Administration (NOAA) Coastal Change Analysis Program (C-CAP) Regional Land Cover layer (NOAA, 2015) to represent most land cover categories, and the TIGER/Line roads layer (USCB, 2014) for added spatial resolution in representing roads. Figure 4-1Error! Reference source not found. Error! Reference source not found. shows NOAA land use coverage (top panel) and NOAA pervious/impervious cover with TIGER roads overlay (bottom panel) for St. Thomas. summarizes the composite land use distribution within TMDL watersheds. Appendix D – Modeling Report June 2022 9 Figure 4-1. NOAA 2015 C-CAP land use/land cover (top panel) and NOAA pervious/impervious cover with TIGER Census roads layer (bottom panel) for Saint Thomas, USVI. Appendix D – Modeling Report 10 June 2022 Table 4-1. Composite land use categories and model groups (Source: 2015 NOAA C-CAP) Model Groups Land Use Category Area (acres) Percent Roads Impervious Roads 1,136.06 4.0% Urban Impervious Impervious 1,571.06 5.6% Urban_Pervious Developed, Open Space 2,173.88 7.7% Agriculture Cultivated Crops 249.29 0.9% Pasture/Hay 1,921.24 6.8% Grass_Shrub Grassland/Herbaceous 718.64 2.6% Forest Deciduous Forest 9,078.06 32.3% Evergreen Forest 1,184.91 4.2% Grass_Shrub Scrub/Shrub 8,748.85 31.1% Wetland_Palustine Palustrine Forested Wetland 273.42 1.0% Palustrine Scrub/Shrub Wetland 2.37 0.0% Palustrine Emergent Wetland 25.47 0.1% Wetland_Estuarine Estuarine Forested Wetland 221.46 0.8% Estuarine Scrub/Shrub Wetland 63.83 0.2% Estuarine Emergent Wetland 5.88 0.0% Barren Unconsolidated Shore 76.17 0.3% Bare Land 313.33 1.1% Water_Open Open Water 341.01 1.2% Total: 28,104.96 100.0% 4.1.2 Soils Data Soils data for the watershed were obtained from the State Soil Geographic and Soil Survey Geographic Database (STATSGO/SSURGO). While the 2005 Summit to Sea Characterization of Coastal Watersheds - US Virgin Islands data set was also identified as a possible source, much of the soil- related information in this dataset was derived from SSURGO data; Therefore, the SSURGO and STATSGO data sets was used directly (NRCS 2016a, NRCS 2016b). The SSURGO data set includes an attributed for hydrologic soil groups (HSG) used to characterize soil infiltration and runoff potential. This attribute is incorporated into the watershed model to parameterize infiltration. HSG is divided into four primary categories using an A through D designation. HSG-A generally has the lowest runoff potential whereas HSG-D has the highest runoff potential. Some map units may contain a compound soil grouping containing multiple HSGs which will be reclassified using a conservative estimate of the lowest infiltrating of the HGSs. The soils database is composed of a GIS layer of polygon map units, and a linked database with multiple tables with soil properties. The representative hydrologic soil group for each polygon (as determined in the STATSGO/SSURGO database) was obtained for each map unit from the map unit aggregated attribute table (MUAGGATT) in the database. Figure 4-2 is a map showing the STATSGO/SSURGO hydrologic soil groups for STEER. Next, the NOAA land use layer was overlaid on the soils layer to further differentiate land use categories by hydrologic soil group. Within a given subwatershed, areas of the same land use category can have different underlying hydrologic soil groups—two areas of the same land use category with different soils could exhibit different runoff Appendix D – Modeling Report June 2022 11 potentials. Table 4-2 summarizes land use area by model group, cross-tabulated with hydrologic soil group distribution. Figure 4-2. STATSGO/SSURGO hydrologic soil groups for St. Thomas, USVI. Table 4-2. Land use and soils overlay (Source Data: 2015 NOAA C-CAP and STATSGO/SURGO Soils) Model Groups Area (acres) Total Hydrologic Soil Group D C B A N/A Roads_Impervious 1,136 4.0% 2.3% 1.5% 0.13% 0.08% 0.00% Urban_Impervious 1,571 5.6% 3.5% 1.8% 0.17% 0.10% 0.00% Urban_Pervious 2,174 7.7% 3.9% 3.3% 0.50% 0.06% 0.00% Agriculture 2,171 7.7% 2.8% 4.5% 0.43% 0.01% 0.00% Grass_Shrub 1 9,467 33.7% 13.7% 18.1% 1.74% 0.13% 0.01% Forest 1 10,263 36.5% 22.8% 11.4% 1.99% 0.31% 0.00% Wetland_Palustine 301 1.1% 0.1% 0.5% 0.14% 0.29% 0.00% Wetland_Estuarine 291 1.0% 0.7% 0.2% 0.01% 0.03% 0.01% Barren 390 1.4% 1.0% 0.3% 0.01% 0.06% 0.03% Water 341 1.2% 0.7% 0.1% 0.00% 0.02% 0.45% Total 28,105 100.0% 51.5% 41.8% 5.12% 1.09% 0.5% 1: Red Box highlights predominant land uses that were further subdivided by Hydrologic Soil Group. Appendix D – Modeling Report 12 June 2022 The overlay summary showed that 93 percent of the study area was composed of either C or D soils. Soil Groups A, B, and unclassified (N/A) made up the other 7 percent. Groups C and D combinations were kept as classified for Forest and Grass_Shrub; however, because the total area was relatively small in comparison to C and D, groups A, B, and N/A were consolidated into a single category to simplify the number of HRU combinations in the model. Soil group is not relevant when overlaid with impervious land use, so all soils for those areas were also consolidated into a single category based on land only. Likewise, the land use categories (Urban_Pervious, Agriculture, Wetland, Barren, and Water), which make up 20 percent of the study area (with most of the area –18 out of the 20 percent of study area—being in soil groups C or D), they were also consolidated by land use only. During model calibration, those areas were considered mixed soils that were parameterized to exhibit an average C/D soils response. 4.1.3 Slope The United States Geological Survey (USGS) National Elevation Dataset (NED) 10-meter digital elevation model (DEM) was used to derive slopes for this analysis. Figure 4-3 shows NED elevation, while Figure 4-4 shows derived slopes. As was done with soils, land use was also overlaid with the DEM-derived slopes to further differentiate certain land use categories by slope where applicable. Table 4-3 is a summary of the land use and slope overlay for the TMDL watersheds. A threshold of 15 percent slope was used to differentiate land between high- and low-slope areas. Areas with slopes <15 percent were classified as “Low”, while areas with slopes ≥15 percent were classified as “High.” Within the TMDL watersheds, 65 percent of the area had “High” slopes, while 35 percent had “Low” slopes. For this study, Grass_Shrub, Forest, and Urban land uses areas, which make up 62 percent of the study area, were further divided into High and Low slope categories. Grass_Shrub and Forest were subdivided because they represent most of the area, and Urban because it represents anthropogenic activity. Subdividing these categories added more resolution for characterizing variability in hydrologic and water quality process responses in the model. Figure 4-3. National Elevation Dataset coverage for St. Thomas, USVI. Appendix D – Modeling Report June 2022 13 Figure 4-4. National Elevation Dataset slope coverage for St. Thomas, USVI. Table 4-3. Land use and slope overlay (Source Data: 2015 NOAA C-CAP and NED DEM slopes) Model Group Area (acres) Area Average Slope <=15% >15% <=15% >15% Roads 1,136 1.73% 2.32% 5.4% 14.6% Urban_Impervious 1 1,571 2.85% 2.74% 5.6% 13.3% Urban_Pervious 1 2,174 5.88% 1.86% 4.8% 12.6% Agriculture 2,171 7.14% 0.59% 3.8% 14.3% Grass_Shrub 1 9,467 8.02% 25.67% 5.6% 16.2% Forest 1 10,263 5.06% 31.46% 5.1% 18.3% Wetland_Palustine 301 1.03% 0.04% 3.1% 13.3% Wetland_Estuarine 291 0.99% 0.03% 2.5 14.9 Barren 390 0.88% 0.51% 3.3 15.2 Water 341 1.03% 0.19% 3.6 11.1 Total 28,105 35% 65% - - 1: Red Box highlights predominant land uses that were further subdivided by slope. Appendix D – Modeling Report 14 June 2022 4.1.4 Other HRU Refinements Two other refinements were made to HRU layer to help characterize unique features that are important to characterizing water quantity and quality in the USVI: (1) rainwater harvesting and (2) construction. There were no specific layers identified that explicitly provided the information necessary; however, estimates could be derived from other supplemental supporting information and literature searches. These refinements are described in this section. 4.1.4.1 Rainwater Harvesting Since the early 1930s, USVI has had mandatory laws requiring private residences and business to build cisterns for rainwater harvesting. The building code reenacted in 1964 and revised in 1996 has a clause requiring all dwellings not connected to municipal water supplies to construct cisterns for capturing, storing, and using rainwater for domestic water supply (Solomon, 2009). Customer service numbers from USVI Water and Power Authority (WAPA) show that 20 percent of residential households and less than 15 percent of businesses have water pipe connections. The other 80 to 85 percent of addresses rely on cistern water as their main source of domestic water (Solomon, 2009). Because the purpose of this study is to support development of implementation plans, a higher- resolution Light Detecting and Ranging (LiDAR) dataset was used to refine assumptions about certain features on the landscape. LiDAR was used to differentiate rooftops impervious area from other urban impervious area. Because most of the rooftop runoff in USVI is used for domestic water supply, it is not available for watershed hydrology. Differentiating rooftop area from other impervious area makes it possible to quantify that volume of water and the role it plays in the water balance. On the basis of WAPA customer numbers, and estimated rooftop area, it was assumed in the model that 85 percent of rooftop water is not available for runoff because it is captured in cisterns for domestic use. 4.1.4.2 Construction Activity Construction activity is typically one of the most impactful anthropogenic factors in a watershed. Figure 4-7 shows where the driveway of a typical construction site meets the main road in St. John. In most cases, it is a temporary activity that occurs for a season; however, while it is occurring, its impact can vary depending the quantity and quality of erosion mitigation measures that have been put in place. Although no spatial layers were identified that explicitly delineated construction activity, a rough estimate can be made by comparing land use changes between two snapshots in time: the 2007 and the 2012 NOAA C-CAT land use layers. Table 4-4 summarizes land use change within the TMDL watersheds. Appendix D – Modeling Report June 2022 15 According to the NOAA C-CAT land use change dataset, about 985 acres in the TMDL watersheds changed from one category to another between the 2007 and 2012 snapshots. Changes in road area was not explicitly represented in this analysis because the TIGER/line roads layer captured one point in time—any new roads would be implicitly reflected among changes to “Urban” land use. Most “changing” area was Agriculture → Grass_Shrub or Agriculture→ Forest, with some Grass_Shrub → Agriculture. Construction area was measured as a change from any undeveloped land use category to Urban. Some area classified as “Barren” in 2007 became “Urban” in 2012, suggesting that at any given point in time, some construction activity may be exhibited as “Barren” in the land use layer. The net change of area from any land use to “Urban” was about 220 acres, which represents a little less than 1 percent of the modeled study area. Changes were spatially tabulated by subwatershed as represented in the change layer; however, the temporal distribution is unknown. Assuming that construction is evenly distributed over the five years (2007-2012), with an average duration of a construction project as approximately one year, up to 20 percent of the changing area at any given time can be characterized as being under construction. Table 4-4. Land use change matrix (Source Data: 2007 and 2015 NOAA C-CAP land use) Changed From (2007) Changed To (2012), acres 1 Land Use, Area (acres) Grass_ Shrub Agriculture Barren Forest Urban Pervious Urban Impervious Water Wetland_ Palustine Wetland_ Palustine Grass_Shrub 348 --- 127 36 69 81 33 2 0 0 Agriculture 319 231 --- 3 72 7 3 2 0 0 Barren 129 46 50 --- 4 17 11 0 0 0 Forest 112 38 17 9 --- 23 24 1 0 0 Urban Pervious 35 14 0 1 7 --- 12 0 0 0 Urban Impervious 19 8 0 2 4 6 --- 0 0 0 Water 21 0 0 0 0 0 0 --- 3 18 Wetland_Palustine 2 0 0 0 0 0 2 0 --- 0 Wetland_Estuarine 1 0 0 0 1 0 0 0 0 --- Total 985 337 194 52 156 134 87 5 3 18 1: Red Box areas that changed to urban from other land uses between 2007 and 2012. Figure 4-7. Example of failed sediment controls at a construction site (St. John, USVI 2016). Appendix D – Modeling Report 16 June 2022 4.1.4.3 Soil Erodibility The advantage of an HRU-based approach is that it retains much of the resolution of spatial variability for model parameterization at the level of the smallest modeling unit (land unit). This minimizes the need to specify diverse combinations of model parameter groups at the subwatershed level. One of the key attributes of the STATSGO/SSURGO soil layer is the K-factor, which is a measure of soil erodibility (detachment and runoff). K-factor was estimated as an area-weighted average value of the top soil layer for each HRU raster pixel. Clay soils, which are resistant to detachment, have low K values (0.05 to 0.15). Likewise, coarse-textured sandy soils that are easily detached, but are not easily mobilized by runoff, also have low K values (0.05 to 0.2). Soils with moderate silt and loam content have moderate K values (0.25-0.4) because they are moderately susceptible to both detachment and runoff. Soils with high silt content are the most erodible of all soils (K > 0.4), because they are easily detached and are associated with high rates of runoff. K-Factor was cross-tabulated for areas having the same hydrologic soil group, as summarized in Table 4-5. The analysis showed that D soils were moderately erodible, while C soils were more erodible, suggesting that C soils should be parameterized as slightly more erodible than D soils, even though D soils produce more runoff than C soils. Soils are dramatically altered by stripping or compacting during construction activities. A survey of native soils at the construction site is usually required for most erosion control plans. Hydrologic soil group (which can also serve as a surrogate for soil erodibility), is estimated based on soil texture, as shown in Table 4-5 and illustrated in Figure 4-5. Within a given soil group, K factor varies with soil texture; however, within the TMDL watersheds, 83 percent of the area is either C or D soils with area- weighted average K-factors of 0.3 and 0.2, respectively. The spatial variability of K-factor is shown in Figure 4-6. Table 4-5. K-Factor vs. hydrologic soil group and associated soil textures (Source Data: NRCS 2016a, 2016b) Hydrologic Soil Group  Low Moderate High → K-Factor (Soil Erodibility) 0 0.1 0.2 0.3 0.4 Area (acres) A 39 0 268 0 0 B 0 0 1,439 0 0 C 0 0 1,488 10,266 0 D 54 237 13,058 922 190 N/A 144 0 0 0 0 Percentage A 0.14% 0.00% 0.95% 0.00% 0.00% B 0.00% 0.00% 5.12% 0.00% 0.00% C 0.00% 0.00% 5.29% 36.53% 0.00% D 0.19% 0.84% 46.46% 3.28% 0.68% N/A 0.51% 0.00% 0.00% 0.00% 0.00% Associated Soil Textures A Sand, Loamy Sand, or Sandy Loam B Silt, Silt Loam or Loam C Sandy Clay Loam D Clay Loam, Silty Clay Loam, Sandy Clay, Silty Clay, or Clay Appendix D – Modeling Report June 2022 17 Figure 4-5. Standard United States Department of Agriculture (USDA) Soil Triangle with Hydrologic Soil Group mapping. Figure 4-6. STATSGO/SSURGO K-Factor for St. Thomas, USVI. Appendix D – Modeling Report 18 June 2022 4.1.4.4 Effective Impervious Area Within a particular modeled subwatershed, HRU segments are modeled as being parallel to one another. Each segment flows directly to the routing segment without any interaction with neighboring segments. On physical landscapes, impervious land will flow downhill over pervious land en route to the watercourse on the landscape. For modeling purposes, Effective Impervious Area (EIA) represents the portion of total, or Mapped Impervious Area (MIA), that is directly-connected to the stream routing segments. Impervious areas that are not connected to the drainage network have the opportunity to flow onto pervious surfaces, infiltrate, and become part of pervious subsurface and overland flow. Because segments are modeled as being parallel to one another in LSPC, this process can be approximated using a conversion of a portion of impervious land to pervious land. On the open landscape, runoff from disconnected impervious surfaces can overwhelm the infiltration capacity of adjacent pervious surfaces during large rainfall/runoff events creating sheet flow over the landscape— therefore, the MIA→EIA translation is not actually a direct linear conversion. Finding the right balance between MIA and EIA can be an important part of the hydrology calibration effort. An analogous parameter from the Cadmus 2011 SUSTAIN runoff modeling study is the fraction of impervious area that is directly-connected to a watercourse (DCIA). The complementary fraction of DCIA (1 – DCIA) is the amount of impervious area within a catchment that is first routed to pervious land before it reaches the watercourse. Among their three selected calibration watersheds, Cadmus assumed that between 80 and 90 percent of mapped impervious area is first routed to pervious area, meaning 10 to 20 percent is DCIA (Cadmus, 2011). During model LSPC calibration for this study, a strong correlation was found between density of developed area and assumed DCIA, which was corroborated by improvements in model performance (reduction in the flashiness of storm peaks). Table 4-6 summarizes total, urban, MIA and EIA, along with calibrated DCIA ranges among the TMDL watersheds. The relationship between impervious density and DCIA was also intuitively reasonable because as developed area density increases in a subwatershed, so does the likelihood of that impervious area somehow being directly connected to a watercourse. Conversely, the less developed a subwatershed is, the potential impact of DCIA also becomes smaller. Table 4-6. Mapped Impervious Area (MIA) vs. Effective Impervious Area (EIA) summarized by TMDL watershed TMDL Watershed Area (acres) Average DCIA2 Range of DCIA (by Subwatershed) Total Total Urban MIA1 EIA3 Min 25th Median 75th Max STEER 3,715 1,182 820 261 32% 3% 17% 33% 40% 71% 1: Mapped Impervious Area (MIA) 2: Directly-Connected Impervious Area (DCIA) Percentage 3: Effective Impervious Area (EIA) = MIA ×DCIA Appendix D – Modeling Report June 2022 19 4.2 Source Characteristics TMDL critical conditions include both temporal and spatial influences. Four primary sources of pollutants were considered: (1) flow and loads from watersheds, (2) point source discharges, (3) septic system discharges, and (4) direct discharges from live-aboard vessels. Previous studies have characterized flow and loads from watersheds. These are shown in Table 4-7. For the temporal component, point sources and septic systems were modeled as constant while watershed and live- aboard vessel discharges varied monthly. For monthly variable sources, long-term continuous simulation loads from the watershed were aggregated and expressed as daily loads per month, while the estimated number of live-aboard vessel residents per month and per-capita boat effluent volume/loads were multiplied to estimate monthly loads from live-aboard vessels. For the spatial component, the critical-condition mixing volume was estimated as a ratio of total watershed flow to equivalent ocean water. This approach implicitly accounts for the influence of physical factors affecting the mixing efficiency of each waterbody. Appendix D – Modeling Report 20 June 2022 Table 4-7. Summary of regional and national unit-area loads by land use (LU) from various literature sources Pollutant Source Land use or Location of Study Selected Model LU for Comparison Minimum Maximum Units ENTERO East End Watersheds Management Plan (2011) Forested Forest 12 12 Billions/acre/yr East End Watersheds Management Plan (2011) Rural Grass 39 39 Billions/acre/yr TSS Ramos-Scharron (2010) Parguera Disturbed Construction 18 26,765 lbs/ac/yr East End Watersheds Management Plan (2011) Forested Forest 50 78 lb/acre/yr Ramos-Scharron (2010) Parguera Undisturbed Forest 1 45 lbs/ac/yr East End Watersheds Management Plan (2011) Rural Grass 100 127 lb/acre/yr TN Reckhow et al. (1980) Non Row Crops Agriculture 0.86 6.96 lb/ac/yr Dodd, R. C., McMahon, G. and Stichter, S. (1992) Developed Developed 4.45 8.65 lb/ac/yr Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Residential Developed 4.45 6.50 lb/ac/yr Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Forest Forest 0.89 5.61 lb/ac/yr Dodd, R. C., McMahon, G. and Stichter, S. (1992) Forest/wetland Forest 0.61 3.38 lb/ac/yr East End Watersheds Management Plan (2011) Forested Forest 2.50 2.51 lb/acre/yr Reckhow et al. (1980) Forested Forest 1.23 5.57 lb/ac/yr Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Idle land Grass 0.45 5.34 lb/ac/yr East End Watersheds Management Plan (2011) Rural Grass 4.60 4.61 lb/acre/yr TP Reckhow et al. (1980) Non Row Crops Agriculture 0.09 2.58 lb/ac/yr Dodd, R. C., McMahon, G. and Stichter, S. (1992) Developed Developed 0.40 1.34 lb/ac/yr Appendix D – Modeling Report June 2022 21 Pollutant Source Land use or Location of Study Selected Model LU for Comparison Minimum Maximum Units Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Residential Developed 0.69 1.96 lb/ac/yr Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Forest Forest 0.01 0.78 lb/ac/yr Dodd, R. C., McMahon, G. and Stichter, S. (1992) Forest/wetland Forest 0.08 0.19 lb/ac/yr East End Watersheds Management Plan (2011) Forested Forest 0.20 0.21 lb/acre/yr Reckhow et al. (1980) Forested Forest 0.02 0.74 lb/ac/yr Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989) Idle land Grass 0.04 0.22 lb/ac/yr East End Watersheds Management Plan (2011) Rural Grass 0.70 0.71 lb/acre/yr BOD5 Mangrove Lagoon and Benner Bay DO TMDL (2003) Agriculture Agriculture 29 29 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Residential (High Density) Developed 26 26 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Residential (Low Density) Developed 14 14 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Residential (Medium Density) Developed 20 20 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Urban Developed 206 206 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Parks/Recreation/ Open Space Grass 6 6 lb/acre/yr Mangrove Lagoon and Benner Bay DO TMDL (2003) Undeveloped Grass 6 6 lb/acre/yr Appendix D – Modeling Report 22 June 2022 4.3 Data Availability and Quality Data from a range of sources were collected and applied for model development. Table 4-8 presents a summary of the quantity and quality of available data that were used for model development. Many of these data sources are described further in context of how they were applied in this report. For example, meteorological (Appendix E), hydrology, and land cover data are key components of the watershed model. Other data sources not directly associated with the watershed are individually described in the context of their respective TMDL documents. Water quality data sheets are summarized in Appendix B and are considered again for the mixing model relative to corresponding watershed model loads. Surrogate pollutants were also identified for some of the TMDL listings, as described in this section. Supporting trends analysis for all listed pollutants are presented in Appendix C. Table 4-8. Data quantity and quality inventory for model development Legend:  Available, Comprehensive  Available, Partial  Available, Deficient - Not Applicable or Available Supporting Information Data Source Data Type or Method Category Information Observed/ Digital Literature/ Surveys Estimated/ Empirical Weather Precipitation NCDC   Temperature NCDC   Evapotranspiration NCDC  Wind Speed NCDC  Hydrology Elevation USGS DEM  Flow records USGS Daily Flow Records  NHD USGS NHD  Stream geometry Estimated   Soils STATSGO/SSURGO  Water Quality Ambient monitoring STORET  Storm monitoring - Groundwater USGS  Land Use/Cover NOAA C-CAP  Hydraulics Tidal information NASA - - - Bathymetry NOAA  Sources Septic numbers/location Census  Point sources NPDES ICIS  MS4 Census Urban Areas  Appendix D - Modeling Report June 2022 23 Legend:  Available, Comprehensive  Available, Partial  Available, Deficient - Not Applicable or Available Supporting Information Data Source Data Type or Method Category Information Observed/ Digital Literature/ Surveys Estimated/ Empirical Live-aboard vessels Surveys   Land Use Change Land Use Layer  Road network Census Tiger Roads  Tourism Rates Bureau of Economic Research   Building Permits Bureau of Economic Research   4.3.1 Meteorological Data Meteorological data are a critical component of the watershed model. They are the forcing functions that drive both hydrology and associated water quality responses. For modeling Caribbean climate conditions, precipitation and evapotranspiration are required at a minimum. In general, hourly (or finer-resolution) precipitation data are recommended for simulating hydrology, erosion, and water quality processes in LSPC. The two primary data sources of observed weather data evaluated and processed for modeling were: (1) the National Climatic Data Center (NCDC) Global Historic Climatology Network (GHCN) daily gage network (30 gages) and (2) the NCDC Quality Controlled Local Climatological Data (QCLCD) hourly gage network (2 gages). Figure 4-7 shows the location of NCDC meteorological stations in the USVI. Figure 4-7. NCDC meteorological gaging stations on St. Thomas, USVI. Appendix D – Modeling Report 24 June 2022 4.3.1.1 Quality Control The Virgin Islands watershed model runs at an hourly time step. Previous experience has shown that the NCDC daily precipitation dataset tends to be more reliable, in terms of total reported volumes, than the hourly dataset; however, hourly gages are still incorporated into the process for rainfall disaggregation. Many of the rainfall stations contain intervals of missing, deleted, or accumulated data. Missing or deleted intervals are periods during which either the gage malfunctioned, or the data records were lost. Accumulated intervals contain cumulative precipitation reported over several hours or days, but the exact temporal distribution of the data is unknown due to a gage malfunction. Two commonly used estimation techniques for patching missing rainfall data are the Normal Ratio Method (Dunne and Leopold, 1978) and the Distance Power Method. The Normal Ratio Method corrects for orographic variability through normalization; however, the Distance Power Method does not consider that. A hybrid approach was used for quality-controlling the Virgin Islands dataset. First, candidate stations were selected from among the gages using nearest distance to gage (a minimum of 3 nearby stations). The search radius was incrementally widened to add stations until at least one good candidate station from the same island was available for each impaired month in the historical record for the station. Stations were preferentially selected, but not weighted, by distance. Second, the Normal Ratio Method was used patch each station with its respective set of nearby stations. Patching of missing and deleted data intervals was performed on a daily time step basis. Next, those patched daily intervals were disaggregated to hourly using hourly rainfall distributions at nearby gages. The Normal Ratio Method estimates missing daily rainfall using a weighted average from surrounding stations with similar rainfall patterns according to the relationship: 𝑃𝐴= 1 𝑛(∑𝑁𝐴 𝑁𝑖 𝑛 𝑖=1 × 𝑃𝑖) where PA is the missing precipitation value at station A, n is the number of surrounding stations with valid data for the same day, NA is the long-term average monthly precipitation at station A, Ni is the long-term average monthly precipitation at nearby station i, and Pi is the observed daily precipitation at nearby station i. For months where data were either zero or heavily impaired (i.e. more than 50% missing), the long-term average annual precipitation value was used for NA and Ni instead of the impaired monthly values. For each missing day at station A, n consists of only the surrounding stations with valid data; therefore, for each day, n can vary from 1 to the maximum number of surrounding stations. When no precipitation is available at the surrounding stations from the same island, zero precipitation is also assumed at station A. In general, gages located at airports and relatively high population/traffic locations tend to have better data quality than those located in more remote locations. Some of the gages contained long time intervals with no rainfall and no missing flags. Sometimes rainfall occurred at nearby stations during the un-flagged missing intervals. Some intervals were as short as one month, while others extended beyond a year. The longer the period of time, the less likely that there would truly be no rainfall occurring during that interval. For this reason, the original record was further adjusted to flag all “dry” intervals greater than 30 days as “missing,” thereby giving them the opportunity to be patched using nearby gages as applicable. Table 4-9 summarizes the final set of 30 rainfall gages from among the three islands considered for modeling. Quality control summaries for each gage are presented in Appendix E. The stations are sorted from highest quality to lowest quality. Attribute information includes elevation, rainfall evaluation period, percent missing (reported and adjusted), and average annual rainfall (reported and adjusted). Appendix D - Modeling Report June 2022 25 Table 4-9. Summary of Virgin Islands NCDC rainfall data compiled and processed for modeling Quality Rank Station Name Station ID Elevation (ft) Evaluation Period1 Percent Missing2 Rainfall (inches/year) First Last Reported Adjusted Reported Adjusted 1 CHRISTIANSTED AP 671624 61 1/2/1980 12/31/2020 0% 2% 36.23 36.94 2 EAST HILL 672560 120 1/5/1980 12/31/2020 0% 7% 38.47 41.12 3 MONTPELLIER 674900 200 2/2/1980 12/31/2020 0% 8% 33.75 36.67 4 CHARLOTTE AMALIE CYRIL 671640 20 3/24/1980 12/31/2020 0% 11% 36.17 39.53 5 BETH UPPER NEW WORKS 670480 110 1/2/1980 12/31/2020 3% 15% 23.80 27.75 6 WINTBERG 679450 645 9/3/1980 12/31/2020 14% 16% 35.71 42.11 7 EAST END 672551 150 1/6/1980 12/31/2020 0% 17% 31.48 36.52 8 CORAL BAY 671790 30 3/4/1981 9/30/2011 0% 18% 27.21 32.71 9 CRUZ BAY 671980 8 1/5/1980 12/31/2020 13% 18% 38.66 46.47 10 REDHOOK BAY ST THOMAS 677600 4 3/3/1980 12/31/2020 0% 24% 18.54 23.36 11 CHRISTIANSTED FORT 671740 30 1/3/1980 12/31/2020 0% 29% 26.19 36.55 12 ESTATE THE SIGHT 672870 130 2/3/1980 6/28/2002 35% 38% 22.16 37.01 13 GRANARD 673677 65 1/5/1980 6/19/2003 38% 38% 23.62 38.34 14 COTTON VALLEY 2 671810 140 2/2/1982 7/31/2001 44% 47% 14.97 28.96 15 ANNALY 670240 700 1/2/1980 2/21/2003 36% 47% 25.61 47.95 16 CANEEL BAY PLANTATION 671316 60 2/3/1980 2/19/1998 50% 53% 20.28 42.35 17 CATHERINEBURG 671348 845 1/5/1980 10/27/1996 54% 57% 19.42 43.28 18 ESTATE FORT MYLNER 672823 200 1/3/1980 8/30/1995 57% 61% 20.89 51.88 19 FOUNTAIN 673150 250 1/2/1980 4/29/1992 67% 67% 15.05 45.30 20 HAM BLUFF LIGHTHOUSE 673880 80 1/23/1980 10/29/1992 65% 68% 8.50 25.57 21 DOROTHEA AES 672440 800 1/4/1980 5/14/1993 64% 68% 13.51 38.50 22 LAMESHUR BAY 674820 170 1/10/1980 3/29/1989 76% 78% 8.46 36.26 23 WATER ISLE 679222 100 1/9/1980 7/25/1990 72% 79% 8.15 38.27 24 FREDERIKSTED 1 SE 673220 80 1/2/1980 8/30/1989 75% 87% 5.72 41.75 25 ESTATE HOPE 672830 390 5/7/1984 8/30/1989 85% 87% 6.70 50.80 26 PROSPECT HILL 676000 214 9/11/2012 12/31/2020 91% 92% 4.57 54.48 27 ESTATE RUST-OP-TWIST 672860 50 3/22/1980 7/29/1984 91% 92% 3.72 41.90 28 BORDEAUX MOUNTAIN 670820 1,110 1/6/1980 7/30/1984 89% 93% 3.14 42.73 29 CHARLOTTE AMALIE 2 671625 15 1/9/1980 11/21/1981 97% 97% 1.51 43.04 30 TAGUE BAY 678621 30 1/7/1980 5/28/1981 98% 99% 0.70 42.24 1. Evaluation period (for data processing): 1/1/1980 – 12/31/2015. Actual start and end date ranges of the gages may extend outside this window. 2. Percent missing: “Reported” are flagged missing intervals in the dataset over the evaluation period. “Adjusted” includes un-flagged dry intervals longer than 30 days. Appendix D – Modeling Report 26 June 2022 Figure 4-8 presents a 25-year summary of annual precipitation totals at Charlotte Amalie Airport (671640) for the period October 1, 1990, through September 30, 2015. The thinner line (dark blue) shows measured precipitation, the thicker line (orange) shows the processed or “patched” precipitation record, and the bar graph (blue-gray) shows percent missing per year. The original dataset had an un-flagged dry interval between 5/9/1996 and 8/6/1998, longer than 2 years. Figure 4-9 shows that Christiansted Airport (671624), which had the best quality data, had a 33-day un-flagged dry interval in 2005 between 4/8 and 5/11. Patching rainfall from nearby gages added about 5 inches of rainfall during that interval, raising the annual total from about 33 inches to 38 inches. Figure 4-8. Total Precipitation at CHARLOTTE AMALIE CYRIL (671640), 1991-2015. Figure 4-9. Total Precipitation at CHRISTIANSTED AP (671624), 1991-2015. Figure 4-10 summarizes patching results from all stations. Annual average rainfall totals confirm just how variable rainfall can be at different locations among the three islands. The top 10 highest quality 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0 10 20 30 40 50 60 70 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 Percent Missing (Bar Graph) Total Precipitation (in/year) 671640: Processed Total Precipitation (in/year) 671640: Measured Total Precipitation (in/year) 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0 10 20 30 40 50 60 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 Percent Missing (Bar Graph) Total Precipitation (in/year) 671624: Processed Total Precipitation (in/year) 671624: Measured Total Precipitation (in/year) Appendix D - Modeling Report June 2022 27 stations are specifically called out in the graph. Those were used for long-term simulation runs for the TMDL; however, during calibration, some of the other gages were selectively assigned for certain periods of time because in some instances (for specific storms), where they had good data, the storms were more representative of localized rainfall in those watersheds. Total estimated rainfall was also plotted against gage elevation to check for any discernable trends. In general, the lowest elevation gages reported the lowest rainfall totals; the totals also somewhat increased with elevation. However, the irregularity in that trend corroborates the belief that other factors like aspect, terrain-induced rain shadow, and island-to-island rainfall variability also influences long-term average rainfall totals. Figure 4-10. Regional NCDC rainfall totals vs. gage elevation. 4.3.2 Potential Evapotranspiration Evaporation and evapotranspiration (ET) account for a large portion of the hydrologic water budget in most places. Simplified methods such as the Hamon and Jensen approaches are often applied because it is convenient to estimate ET as a function of fewer meteorological constituents. Because of the relatively narrow seasonal variability in climate conditions, those methods are often as robust as the more rigorous energy-balance based approaches. The Pathogen Study used the Hamon method with default coefficient values to derive ET for SUSTAIN (Cadmus, 2011). For this effort, ET was also computed using the default Hamon method, which estimates ET as a function of temperature. Daily ET is disaggregated to hourly over daylight hours as a function of latitude. Potential ET varies as a function of vegetative cover. ET coefficients were applied as a function of land cover to reflect stratification of vegetative impacts, as summarized in Table 4-10. This approach provides more spatial resolution by allowing it to vary as a function of land use, which varies by subwatershed. Two gages among the islands had the best quality long-term temperature timeseries for computing ET: (1) Charlotte Amalie Cyril E. King Airport on St. Thomas, and (2) Christiansted Hamilton Field Airport on St. Croix. The Cyril E. King Airport was used for the northern islands of St. Thomas and St. John, while Hamilton Field Airport was used for St. Croix. Although each island is assigned one timeseries as the basis for ET, land-cover-specific stratification (Table 4-10) adds spatial variability as a function of vegetative cover. Appendix D – Modeling Report 28 June 2022 Table 4-10. Estimated stratification of modeled ET by land cover type Cover Type Land Cover Evapotranspiration Multiplier Rationale 1 Urban Impervious 1.2 Above average ET (warm exposed surfaces) Pervious 0.9 Grass or shrub vegetation Construction 1.0 No vegetation, use standard ET rate Rural Agriculture 0.9 Grass or shrub vegetation Barren 0.9 Grass or shrub vegetation Forest/Wetland 0.85 Light wind, high relative humidity Grass-Shrub 0.9 Grass or shrub vegetation Water Water 1.0 Use evaporation rate for open water 1. Reference: Bedient and Huber, 2002. Table 1.2, Page 47. 4.4 Data Gaps 4.4.1 Hydrological Data Gaps Aside from the meteorological data gaps described in Section 4.3.1, there were very few USVI flow gages with high-quality long-term continuous flow monitoring. For the reported gage data, there was also a relatively high frequency of “estimated” flags among the USGS flow data, suggesting that the ephemeral-to-intermittent flow conditions in the stream channels (commonly referred to as “guts”) have been challenging to measure accurately. In fact, USGS flow monitoring in USVI ended altogether in 2006 due to government budget shortfalls. Table 4-11 summarizes the quantity and quality of USGS data in USVI. The USGS StreamStats application was used to develop flow statistics for the selected calibration sites to evaluate model performance during the more recent period from 2006 through 2020 when observed, continuous data was not available. 4.4.2 Water Quality Data Gaps There were no instream monitoring stations within the guts of any of the watersheds. All of the gages were in the estuarine portions of listed TMDL waterbodies. This data gap was addressed by setting water quality parameter values in the watershed model based on literature values shown in Table 4-7. Appendix D - Modeling Report June 2022 29 Table 4-11. Quantity and quality of USGS streamflow records by water year (Source Data: USGS) Data Quantity (Blue: Percent Complete) and Data Quality (Red: Percent Estimated) Legend (% Complete/Estimated): ○ 0% ◔ 10%-20% ◑ 30%-60% ◕ 70%-90% ● 100% USGS Station ID Water Years (October 1, 1985 – September 30, 2003) 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 50274000 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ◔ ○ ○ ○ ◔ ○ ◔ ◔ ◑ ◕ ◔ ◔ ○ ◔ 50292600 ● ● ○ ○ 50294000 ◔ ● ● ○ ○ ◔ 50295000 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ○ ○ ○ ○ ○ ○ ○ ◔ ◔ ○ ◔ ○ ◔ ◔ ◑ ◔ ◔ ◑ 50295500 ◔ ◕ ○ ○ 50332000 ◕ ● ○ ○ 50333500 ● ● ● ◑ ○ ○ ○ ○ 50334500 ◕ ◕ ◕ ● ● ◔ ○ ○ ○ ○ ○ ○ 50337500 ◑ ◑ ○ ○ 50345000 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ◑ ◔ ◔ ◔ ○ ○ ○ ○ ○ ◔ ◔ ◑ ◔ ◔ ◕ ◔ ◑ ○ ○ ○ ○ 50348000 ◑ ◑ ○ ○ 50349000 ◕ ◔ ○ ○ Appendix D – Modeling Report 30 June 2022 4.5 Time Frame of Simulation The model simulation and climate inputs cover the period from October 1, 1980 through September 30, 2020. Calibration was conducted from 10/1/1992 through 9/30/2006 to coincide with the most recent continuous streamflow data. Additional validation was performed using available streamflow gages statistics for a more recent period from 10/1/2011 through 9/30/2020. 5 MODEL CALIBRATION Calibration refers to the configuration and refinement of the analytical instruments that are used to generate analytical data. As alluded to in the Quality Assurance Project Plan (QAPP) document for this project, by analogy the “instrument” is the predictive tool (i.e., the model) that is to be developed and/or applied. The goal is to link delineated characteristics associated with model segmentation with observed conditions the water bodies of interest. Figure 5-1 is a schematic describing a process for model calibration that aims to minimize the propagation of uncertainty. This process builds upon the model development cycle and elements of data quality control previously described in Section 2.1. Figure 5-1. Process for model calibration to minimize propagation of uncertainty. This study builds upon the efforts and lessons learned from previous modeling studies in the region. As described by Cadmus (2011): Previous studies of USVI hydrology have identified minimal overland flow in vegetated upland regions, and point to rapid lateral subsurface flow and saturation overland flow in low-lying areas as primary runoff sources (Macdonald 1997; Ramos-Scherron & Macdonald 2007). Delayed groundwater outflow (baseflow) is highly variable by location and, where present, occurs during a short portion of the year, resulting in ephemeral to intermittent flow conditions in stream channels (U.S. Geological Survey, 1996). Appendix D - Modeling Report June 2022 31 The model calibration process aims to establish causal relationships between meteorological conditions, delineated model components, and predicted hydrology and water quality responses, in comparison to observed conditions. For this reason, using data for time intervals at locations with the best observed data quality minimizes the chance of propagating error associated with data uncertainty. In keeping with a top-down weight-of-evidence-based methodology, model calibration focused on specific time intervals at locations with overlapping periods of high rainfall and streamflow data quality. At each location, parameters associated with the predominant drainage area HRU characteristics (i.e., land use, slope, or soil group) were the focus of calibration. Table 5-1 summarizes upstream drainage area HRU characteristics for each of the USGS gages presented in Table 4-11. In general, smaller upstream gages, which were predominantly high-sloped forested watersheds, were first used to calibrate the relative magnitude and duration of background flow conditions. The larger downstream watersheds (which exhibited mixed land use, soils, and slope characteristics) where then used to adjust other HRU parameters. Table 5-1. Upstream drainage area HRU characteristic for each USGS gage USGS Station ID Model Drainage Area (acres) Land Use Slope Hydrologic Soil Group Forest Wetland Grass Shrub Agriculture Barren Developed High (≥15%) Low (< 15%) A/B C D Mixed 50274000 1,490 46% 20% 1% 33% 36% 64% 3% 0% 62% 35% 50292600 249 99% 0% 0% 1% 90% 10% 0% 91% 6% 2% 50294000 947 93% 1% 0% 6% 56% 44% 2% 2% 89% 7% 50295000 240 81% 3% 1% 16% 60% 40% 0% 0% 84% 16% 50295500 118 75% 5% 0% 20% 56% 44% 0% 0% 79% 21% 50332000 922 52% 26% 3% 18% 34% 66% 48% 4% 27% 20% 50333500 2,526 36% 21% 20% 23% 11% 89% 4% 16% 36% 44% 50334500 2,666 31% 28% 20% 20% 16% 84% 15% 25% 19% 41% 50337500 137 66% 34% 0% 0% 73% 27% 0% 0% 99% 1% 50345000 1,372 68% 26% 2% 4% 50% 50% 0% 3% 92% 6% 50348000 231 56% 36% 0% 8% 64% 36% 0% 65% 27% 8% 50349000 59 15% 76% 0% 9% 80% 20% 0% 82% 9% 9% 5.1 Water Balance Calibration Calibration of the modeled water balance primarily involved the first two elements in Figure 5-1 (Weather Data and Land Hydrology). Key model components influencing hydrology, hydraulics, and the overall water balance that were evaluated as part of model configuration included: (1) precipitation data quantity and quality, (2) evaporation and evapotranspiration rates. Rainfall in USVI is often localized—two nearby gages can often report notably different rainfall totals for the same storm event. When results did not match observed trends, wind data and historical hurricane logs were referenced to qualify the direction and movement of storms. Rainfall gages were initially assigned to modeled subwatersheds using the Thiessen polygon method. This approach is typically applied to stations that are several miles apart and cover a wide watershed area; however, because an entire island is sometimes 6 to 15 miles from one side to another, any gage Appendix D – Modeling Report 32 June 2022 on the island could be considered as representative depending on the movement of specific storms. Because the islands are small and storms can be variable and localized, the default station assigned using the Thiessen polygon method did not always represent the storms for all time periods as well as different gage. For this reason, different gages on each island were tried at each calibration point to see which one was the most representative of most storms. Performance varied from storm to storm— some were better reflected at certain gages than at others. For each calibration watershed, the rainfall gage on each island that captured most of the storms was ultimately selected as the gage to be used for model calibration. The top-down weight of evidence approach was implemented in the following sequence (1) calibrate background conditions, (2) add intermediate mixed land use areas, and (3) aggregate all sources via routing to a downstream location for comparison with co-located flow data. Figure 5-2 is a schematic showing the parameterization and calibration sequence for land hydrology. Figure 5-2. Model parameterization and calibration sequence for land hydrology. For hydrologic calibration of HSPF, a variety of performance targets have been specified, including Lumb et al. (1994) and Donigian (2000). As previously noted, the LSPC model is functionally identical to USEPA’s HSPF model in terms of hydrologic and water quality processes. Based on those literature sources, performance targets for simulation of the water balance components are summarized in Table 5-2. The general acceptance criterion for models to be applied in this project is to achieve a quality of fit of “Good” or better. Appendix D - Modeling Report June 2022 33 Table 5-2. Typical performance targets for LSPC hydrology simulation (modeled vs. observed) Model Component Very Good Good Fair Poor Error in Total Volume < 5% 5-10% 10-15% > 15% Error in 10% Highest Flow Volumes < 10% 10-15% 15-25% > 25% Calibration outcomes are explained, with some insights and rationale provided for both “Very Good” to “Poor” calibration metrics in light of the top-down weight-of-evidence approach described above. The variety of graphical and statistical points of comparison include: • Time series plots of modeled versus observed streamflow (Section 6.1 and Section 6.2) • Numerical calibration metrics vs. performance targets (Section 6.1, Section 6.2, and Section 6.3) 5.2 Model Parameters Using the approach outlined in the previous section, the final calibrated parameters for LSPC watershed model are presented in Table 5-3. Table 5-3. Summary of calibrated hydrology parameters Parameter Description Unit Value LZSN Lower Zone Nominal Soil Moisture Storage inches 0-0.7 INFILT Index to Infiltration Capacity in/hr 0-0.2 KVARY Variable groundwater recession 1/inches 0-3.0 AGWRC Base groundwater recession none 0-0.87 CEPSC Interception storage capacity inches 0.05-0.15 UZSN Upper zone nominal soil moisture storage inches 0-1.1 NSUR Manning’s n (roughness) for overland flow none 0.05-0.3 INTFW Interflow inflow parameter none 0 IRC Interflow recession parameter none 0-0.8 LZETP Lower zone ET parameter none 0-0.7 PETMAX Temp below which ET is reduced deg. F 45 PETMIN Temp below which ET is set to zero deg. F 35 INFEXP Exponent in infiltration equation none 2 INFILD Ratio of max/mean infiltration capacities none 2 DEEPFR Fraction of GW inflow to deep recharge none 0-0.1 BASETP Fraction of remaining ET from baseflow none 0-0.1 AGWETP Fraction of remaining ET from active GW none 0-0.7 Appendix D – Modeling Report 34 June 2022 5.3 Water Quality Calibration Other factors, besides pollutant loading, influence how a waterbody responds to loadings from various pollutant sources. As previously described, the mixing model methodology uses a mixing ratio as a basis for normalizing the influence of these competing factors. The mixing ratio measures how daily watershed flow volumes of background ocean water mix flow and loads from contributing sources. Geometric analysis of each assessment unit was performed to provide a physical basis for selected mixing ratios. Bathymetric data for each assessment unit were used to calculate (1) average depth, (2) the ratio of bottom depth at the opening to average depth (which indicates efficiency of lateral mixing between the waterbody and open waters) and (3) the ratio of waterbody opening to shoreline perimeter (which indicates how narrow or wide the opening of the waterbody is for mixing). For each of these indicators, the higher the value, the more conducive it is for improving mixing efficiency. For illustrative purposes, a relative symbolic scale was developed for each of the three indicators to illustrate the impact of that indicator on mixing efficiency. Table 5-4 presents a legend of numeric ranges corresponding to the symbols show Table 5-5, which summarizes those geometric properties for each assessment unit. For each waterbody, the mixing ratio was adjusted and calibrated so that diluted model concentrations resembled observed concentrations. The dilution volume was validated for each waterbody relative to others in reference to physical attributes such as width/narrowness of the estuary opening, bathymetric obstructions, orientation of the estuary opening relative to upwelling and wind-driven storm movement, and the location of monitoring gages within the estuary. Critical-condition concentrations were calculated for each pollutant for each representative day per month. The limiting month for each pollutant was the month with the highest modeled pollutant concentration. This result was also used to validate the model. Table 5-4. Legend for geometric properties table Legend of Geometric Properties (See Table 5-5) Entry, Interpretation Average Depth Lateral Mixing Opening Width ● High Over 10 m No bottom obstruction Over 45%  Medium 3 - 10 m Some obstruction 25 - 45% ○ Low 3 m or less Obstructed 25% or less Appendix D - Modeling Report June 2022 35 Table 5-5. Summary of measured geometric properties and calibrated critical mixing ratio and volume for each Assessment Unit Watershed Waterbody ID Waterbody Name Mixing Ratio* Dilution Volume Geometric Properties m3 % of Total Average Depth Lateral Mixing Opening Width Opening Direction St. Thomas East End Reserve (STEER) STT-25 Great Bay 50.0 29,297 0.2% ● ●  N/E STT-28 Cowpet Bay 7.0 3,960 0.5%  ●  S STT-31 Nazareth Bay 8.0 4,526 0.2%  ●  S STT-33 Benner Bay 10.0 8,085 0.3%  ●  S STT-34 Benner Bay Lagoon Marina 8.0 4,526 9.83% ○ ○ ○ S STT-35 Mangrove Lagoon 2.5 9,486 1.9% ○ ○ ○ S Symbol Legend: ● High,  Medium, ○ Low (See Table 5-4 for corresponding numeric ranges for each symbol) * Mixing Ratio: Number of average daily watershed flow volumes used for dilution. Appendix D – Modeling Report 36 June 2022 6 MODEL EVALUATION AND ACCEPTANCE Model evaluation focused on confirming the model’s ability to represent systematic long-term behavior seen in the observed data. This process focused on time periods with overlapping intervals of both high-quality observed precipitation and flow. Due to the limitations in data available across the USVI, gages on islands other than St. Thomas were needed to demonstrate model performance. 6.1 Calibration Example: Jolly Hill Gut, St. Croix At Jolly Hill Gut, daily-intersected rainfall and flow intervals were summarized over a 10-year period (10/1/1996 – 9/30/2006). Table 6-1 summarizes modeled vs. observed results during intervals represented by the four combinations of good and estimated rainfall vs. flow data quality. Model calibration focused on intervals with both “good” observed rainfall and “good” observed flow, as represented by the red-outlined boxes in Table 6-1, which represented 68 percent of both observed flow days and observed flow volume. Figure 6-1 is a monthly summary of modeled vs. observed flow and data quality. The top hyetograph also includes percent missing rainfall (gray), which was estimated from nearby gages using the Normal Ratio Method. The bottom flow time series graph shows both percent estimated (light blue) and percent missing (gray) for observed flow as reported by USGS at the gage. On the basis of this review, a two-year continuous period (10/1/1997 – 9/27/1999) was selected as the most representative calibration time period at Jolly Hill Gut. It included a dry year and a wet year (including the remnants of Tropical Storm Nicole, Nov-Dec 1998). It also represented the longest stretch of continuous “good” observed rainfall and “good” observed flow data. Figure 6-2 through Figure 6-5 show daily, monthly, seasonal average, and seasonal interquartile comparisons of modeled versus observed hydrology for the selected calibration period. Table 6-1. Model performance at Jolly Hill Gut for good-and-estimated rainfall vs. flow data quality intervals Rainfall Data Quality Percent of Days Percent of Flows Percent Difference Flow Data Quality Flow Data Quality Flow Data Quality Observed Modeled Observed Modeled Observed Modeled Observed 68% 23% 68% 24% 7% 17% Modeled 5% 4% 5% 3% -45% 35% Appendix D - Modeling Report June 2022 37 Figure 6-1. Temporal data quality summary for model calibration at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). Figure 6-2. Daily modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0.0 2.0 4.0 6.0 8.0 10.0 12.0 10/1/1996 2/1/1997 6/1/1997 10/1/1997 2/1/1998 6/1/1998 10/1/1998 2/1/1999 6/1/1999 10/1/1999 2/1/2000 6/1/2000 10/1/2000 2/1/2001 6/1/2001 10/1/2001 2/1/2002 6/1/2002 10/1/2002 2/1/2003 6/1/2003 10/1/2003 2/1/2004 6/1/2004 10/1/2004 2/1/2005 6/1/2005 10/1/2005 2/1/2006 6/1/2006 Percent of Flow (Estimated or Missing) Normalized Streamflow (in.) Percent Estimated Percent Missing Observed: JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI Modeled Streamflow 0% 25% 50% 75% 100% 0 10 20 30 Percent Missing Rainfall (inches) Rainfall (inches/month) Percent Missing: EAST HILL (672560) Appendix D – Modeling Report 38 June 2022 Figure 6-3. Monthly modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). Figure 6-4. Seasonal average modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). 0 10 20 30 40 50 60 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 10/1/1997 1/1/1998 4/1/1998 7/1/1998 10/1/1998 1/1/1999 4/1/1999 7/1/1999 Precipitation (in.) Normalized Streamflow (in.) Precipitation Observed: JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI Modeled Streamflow 0 10 20 30 40 50 60 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Precipitation (in.) Normalized Streamflow (in.) Aggregated Monthly (10/01/1997 - 09/30/1999) Precipitation Observed Modeled y = 0.996x R² = 0.969 0 0.5 1 1.5 0 0.5 1 1.5 Model Streamflow (in.) Observed Streamflow (in.) Normalized Monthly Streamflow Appendix D - Modeling Report June 2022 39 Figure 6-5. Seasonal interquartile modeled versus observed streamflow at JOLLY HILL GUT AT JOLLY HILL, ST. CROIX USVI (Station ID: 50345000). 6.2 Calibration Example: Turpentine Run, St. Thomas At Turpentine Run, daily-intersected rainfall and flow intervals were summarized over a 14-year period (10/1/1992 – 9/30/2006). Model calibration focused on time periods with overlapping intervals of both high-quality observed precipitation and flow. Table 6-2 summarizes modeled vs. observed results during intervals represented by the four combinations of good and estimated rainfall vs. flow data quality. Model calibration focused on intervals with both “good” observed rainfall and “good” observed flow, as represented by the red-outlined boxes in Table 6-2, which represented more than 50 percent of both observed flow days and observed flow volume. Within the observed flow record, there was an anomalous 2-day event (9/15-9/16/1995), which neither appeared in the rainfall record as impaired for flow nor rainfall. There was also no reported tropical storm for those days. When those two days were removed from the observed flow record, the percent difference between modeled and observed flow improved from -19 percent to +7 percent. Figure 6-6 is a monthly summary of modeled vs. observed flow and data quality (the anomalous observed flow value is also shown in the graph). Figure 6-7 through Figure 6-10 show daily, monthly, seasonal average, and seasonal interquartile comparisons of modeled versus observed hydrology for the selected calibration period. y = 0.960x R² = 0.975 0 1 2 3 0 1 2 3 Modeled Streamflow (in.) Observed Streamflow (in.) Normalized Monthly Streamflow 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Streamflow (cfs) Aggregated Monthly (10/01/1997 - 09/30/1999) Interquartile Ranges Observed Modeled Appendix D – Modeling Report 40 June 2022 Table 6-2. Model performance for Turpentine Run at Mt. Zion, St. Thomas, for good-and-estimated rainfall vs. flow data quality intervals Evaluation Period Rainfall Data Quality Percent of Days Percent of Flows Percent Difference Flow Data Quality Flow Data Quality Flow Data Quality Observed Modeled Observed Estimated Observed Estimated All Days Observed 53% 26% 61% 26% -19% -46% Modeled 18% 3% 11% 1% -65% -82% Excluding 9/15 and 9/16/1995 Observed 53% 26% 54% 32% 7% -46% Modeled 18% 3% 11% 2% -54% -82% Figure 6-6. Temporal data quality summary for model calibration at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0 10/1/1992 3/1/1993 8/1/1993 1/1/1994 6/1/1994 11/1/1994 4/1/1995 9/1/1995 2/1/1996 7/1/1996 12/1/1996 5/1/1997 10/1/1997 3/1/1998 8/1/1998 1/1/1999 6/1/1999 11/1/1999 4/1/2000 9/1/2000 2/1/2001 7/1/2001 12/1/2001 5/1/2002 10/1/2002 3/1/2003 8/1/2003 1/1/2004 6/1/2004 11/1/2004 4/1/2005 9/1/2005 2/1/2006 7/1/2006 Percent of Flow (Estimated or Missing) Normalized Streamflow (in.) Percent Estimated Percent Missing Observed: TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI Modeled Streamflow 0% 25% 50% 75% 100% 0 5 10 15 20 Percent Missing Rainfall (inches) Rainfall (inches/month) Percent Missing: CHARLOTTE AMALIE CYRIL (671640) Appendix D - Modeling Report June 2022 41 Figure 6-7. Daily modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). Figure 6-8. Monthly modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). 0 2 4 6 8 10 12 14 16 18 20 0 20 40 60 80 100 120 140 10/1/1996 2/1/1997 6/1/1997 10/1/1997 2/1/1998 6/1/1998 10/1/1998 2/1/1999 6/1/1999 10/1/1999 2/1/2000 6/1/2000 10/1/2000 2/1/2001 6/1/2001 10/1/2001 2/1/2002 6/1/2002 10/1/2002 2/1/2003 6/1/2003 10/1/2003 2/1/2004 6/1/2004 10/1/2004 2/1/2005 6/1/2005 10/1/2005 2/1/2006 6/1/2006 Precipitation (in.) Streamflow (cfs) Precipitation (in.) Observed: TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI Modeled Streamflow 0 10 20 30 40 50 60 0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 10/1/1996 1/1/1997 4/1/1997 7/1/1997 10/1/1997 1/1/1998 4/1/1998 7/1/1998 10/1/1998 1/1/1999 4/1/1999 7/1/1999 10/1/1999 1/1/2000 4/1/2000 7/1/2000 10/1/2000 1/1/2001 4/1/2001 7/1/2001 10/1/2001 1/1/2002 4/1/2002 7/1/2002 10/1/2002 1/1/2003 4/1/2003 7/1/2003 10/1/2003 1/1/2004 4/1/2004 7/1/2004 10/1/2004 1/1/2005 4/1/2005 7/1/2005 10/1/2005 1/1/2006 4/1/2006 7/1/2006 Precipitation (in.) Normalized Streamflow (in.) Precipitation Observed: TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI Modeled Streamflow Appendix D – Modeling Report 42 June 2022 Figure 6-9. Seasonal average modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). Figure 6-10. Seasonal interquartile modeled versus observed streamflow at TURPENTINE RUN AT MT. ZION, ST. THOMAS USVI (Station ID: 50274000). 6.3 Validation Example: Turpentine Run, St. Thomas An independent validation period from October 1, 2011 through September 30, 2020 was evaluated for the Turpentine Run at Mt. Zion, St. Thomas streamflow gage comparing mean and 95th percentile statistics from the model to comparable values derived using the USGS StreamStats application (USGS, 2021). The mean annual flow statistic was chosen to evaluate the overall water balance represented in the model while the 95th percentile statistics was chosen to specially evaluate wet weather when non- point source runoff and pollutant loading is highest. Table 6-3 presents a comparison of the mean annual and 95th percentile flow between StreamStats and the LSPC watershed model for the Turpentine 0 10 20 30 40 50 60 0.0 0.5 1.0 1.5 2.0 2.5 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Precipitation (in.) Normalized Streamflow (in.) Aggregated Monthly (10/01/1996 - 09/30/2006) Precipitation Observed Modeled y = 0.911x R² = 0.965 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 Model Streamflow (in.) Observed Streamflow (in.) Normalized Monthly Streamflow y = 0.718x R² = 0.746 0 5 10 0 5 10 Modeled Streamflow (in.) Observed Streamflow (in.) Normalized Monthly Streamflow 0.0 0.5 1.0 1.5 2.0 2.5 Oct Nov Dec Jan Feb Mar Apr May Jun Jul Aug Sep Streamflow (cfs) Aggregated Monthly (10/01/1996 - 09/30/2006) Interquartile Ranges Observed Modeled Appendix D - Modeling Report June 2022 43 Run at Mt. Zion, St. Thomas gage. Generally, the model predictions align with the expected statistics and meet the criteria for calibration presented in Table 5-2. Table 6-3. Comparison of the mean annual and 95th percentile flow between StreamStats and model for Turpentine Run at Mt. Zion, St. Thomas Evaluation Period Metric Mean Annual Flow 95th Percentile Flow StreamStats Modeled StreamStats Modeled 10/1/2011 – 9/30/2020 Flow (cfs) 2.9 2.5 3.9 3.5 % Error -- -13.8% -- -10.3% 6.4 Water Quality Model Observed data consistently showed that November has the highest bacteria concentrations in most of the water bodies. November also corresponds with the beginning of the season when the live-aboard population increases in USVI. Table 6-4 presents the critical condition month with the highest concentration for each modeled pollutant. Table 6-4. Critical condition month for each TMDL waterbody and pollutant under existing conditions Watershed Waterbody ID Waterbody Name Critical-Condition Month TSS BOD5 TN TP ENTERO St. Thomas East End Reserve (STEER) VI-STT-34 Benner Bay Lagoon Marina 10 1 2 2 10 VI-STT-35 Mangrove Lagoon 8 7 7 7 9 VI-STT-25 Great Bay 10 7 10 10 10 VI-STT-28 Cowpet Bay 10 3 3 3 10 VI-STT-31 Nazareth Bay 10 3 2 3 10 VI-STT-33 Benner Bay 10 11 10 10 11 Due to the lack of instream water quality monitoring (all gages were in estuarine portions of the listed TMDL waterbodies), stratifying loads across contributing sources in a representative manner was an important modeling objective. First, event-mean concentrations (EMC) for each pollutant were stratified and assigned to each HRU in the model. Second, modeled runoff and pollutant loads were summarized by land use and validated against unit-area loading rates from literature. In the absence of instream monitored water quality from the watershed, the EMC assignments with comparison of locally estimated unit-area loads formed the basis for representing watershed-originated pollutant loads for this study. Some literature estimates were derived from local TMDLs and watershed management plans, but some others were referenced from other studies. Figure 6-11 shows modeled versus literature values for unit-area sediment loading by land use, slope, and hydrologic soil group among all modeled USVI subwatersheds. The blue quartile box and whiskers represent the range of annual unit-area loads from all modeled subwatersheds, while the red boxes outline the range of unit- area loads (minimum to maximum) from among the reviewed literature. Table 4-7 previously summarized the range of unit-area loads from the respective data sources for each pollutant and land use type—the red boxes in the box plots represent the range of these literature values. Modeled vs. literature values for the other pollutants BOD5 (Figure 6-12), ENTERO bacteria (Figure 6-13), TN (Figure 6-14), and TP (Figure 6-15), are subsequently presented. Appendix D – Modeling Report 44 June 2022 Figure 6-11. Modeled unit-area sediment (lbs/acre) for all watersheds. Figure 6-12. Modeled unit-area BOD5 (lbs/acre) for all watersheds. 73 15 36 527 56 748 48 40 10 20 30 165 56 1 10 100 1,000 10,000 100,000 Agriculture, Barren Forest, Wetland Grass, Shrub Roads Pervious & Impervious Construction High (>15%) Low (≤15%) A & B C D Mixed Impervious Undeveloped Land Use Developed Land Use Slope Hydrologic Soil Group Median/Quartiles Literature Values Average Sediment (lbs/acre) 2 1 2 23 5 3 2 2 0 1 2 2 5 0 50 100 150 200 250 Agriculture, Barren Forest, Wetland Grass, Shrub Roads Pervious & Impervious Construction High (>15%) Low (≤15%) A & B C D Mixed Impervious Undeveloped Land Use Developed Land Use Slope Hydrologic Soil Group Median/Quartiles Literature Values Average BOD5 (lbs/acre) Appendix D - Modeling Report June 2022 45 Figure 6-13. Modeled unit-area ENTERO (Billion/acre) for all watersheds. Figure 6-14. Modeled unit-area TN (lbs/acre) for all watersheds. 3 1 2 124 24 4 3 8 1 1 2 2 24 0 20 40 60 80 100 120 140 160 Agriculture, Barren Forest, Wetland Grass, Shrub Roads Pervious & Impervious Construction High (>15%) Low (≤15%) A & B C D Mixed Impervious Undeveloped Land Use Developed Land Use Slope Hydrologic Soil Group Median/Quartiles Literature Values Average Enterococcus (Billion/acre) 2 1 1 8 4 2 2 2 0 1 1 1 4 0 1 2 3 4 5 6 7 8 9 10 Agriculture, Barren Forest, Wetland Grass, Shrub Roads Pervious & Impervious Construction High (>15%) Low (≤15%) A & B C D Mixed Impervious Undeveloped Land Use Developed Land Use Slope Hydrologic Soil Group Median/Quartiles Literature Values Average Total Nitrogen (lbs/acre) Appendix D – Modeling Report 46 June 2022 Figure 6-15. Modeled unit-area TP (lbs/acre) for all watersheds. 0.2 0.0 0.0 0.4 0.2 0.1 0.1 0.1 0.0 0.0 0.0 0.1 0.2 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Agriculture, Barren Forest, Wetland Grass, Shrub Roads Pervious & Impervious Construction High (>15%) Low (≤15%) A & B C D Mixed Impervious Undeveloped Land Use Developed Land Use Slope Hydrologic Soil Group Median/Quartiles Literature Values Average Total Phosphorus (lbs/acre) Appendix D - Modeling Report June 2022 47 7 REFERENCES Bedient, P. and W. Huber. 2002. Hydrology and Floodplain Analysis. Prentice-Hall, Third Edition. Upper Saddle River, NJ. Cadmus. 2011. Watershed Characterization and Planning for Pathogen Source Reduction in the U.S. Virgin Islands. Prepared for U.S. Environmental Protection Agency Region 2. September 2011. DPNR (US Virgin Islands, Department of Planning and Natural Resources, Division of Environmental Protection). 2004. Dissolved Oxygen Total Maximum Daily Load, Salt River Bay, St. Croix, USVI. Prepared by Tetra Tech, Inc. September 21, 2004. Dodd, R. C., McMahon, G. and S. Stichter. 1992. Watershed Planning in the Albemarle-Palmico Estuarine System. Albemarle-Pamlico Estuarine Study, NC Dept. of Environment, Health, and Natural Resources. Donigian, A.S. Jr. 2000. HSPF Training Workshop Handbook and CD. Lecture #19. Calibration and Verification Issues. Prepared for U.S. Environmental Protection Agency, Office of Water, Office of Science and Technology, Washington, DC. Dunne, T. and L. B. Leopold. 1978. Water in Environmental Planning. W.H. Freeman and Company. San Francisco, California. 818 pp. Horsley Witten Group, Inc. St. Croix East End Watersheds Management Plan, US Virgin Islands. November 2011. Loehr, R. C., Ryding, S. O., and Sonzogni, W. C. (1989). Estimating the Nutrient Load to a Waterbody. "The Control of Eutrophication of Lakes and Reservoirs, Volume I, Chapter 7." Lumb, A.M., R.B. McCammon, and J.L. Kittle Jr. 1994. User’s Manual for an Expert System (HSPEXP) for Calibration of the Hydrological Simulation Program–FORTRAN. Water- Resources Investigations Report 94-4168. U.S. Geological Survey, Reston, VA. Ramos-Scharron, Carlos. 2010. Sediment production from unpaved roads in a sub-tropical dry setting - Southwestern Puerto Rico. Catena 82 (2010) 146-158. Reckhow et al. 1980. Modeling Phosphorus Loading and Lake Response under Uncertainty: A Manual and Compilation of Export Coefficients. US Environmental Protection Agency, Washington, DC. NOAA (National Oceanic and Atmospheric Administration) Coastal Services Center. 2015. The Coastal Change Analysis Program (C-CAP) Regional Land Cover. Charleston, SC NOAA Coastal Services Center. Accessed at www.csc.noaa.gov/digitalcoast/data/ccapregional/index.html NRCS (Natural Resources Conservation Service, United States Department of Agriculture). 2016a. Soil Survey Geographic (SSURGO) Database. Available online at https://sdmdataaccess.sc.egov.usda.gov. Accessed 01/29/2016. Appendix D – Modeling Report 48 June 2022 NRCS (Natural Resources Conservation Service, United States Department of Agriculture). 2016b. U.S. General Soil Map (STATSGO2). Available online at https://sdmdataaccess.sc.egov.usda.gov. Accessed 01/29/2016. Solomon, H. 2009. Effects of Rainwater Harvesting and Desalination on Supply and Demand of Public Water in U.S. Virgin Islands. Proceedings: 4th IASME / WSEAS Int. Conf. on Wat. Res., Hydraulics & Hydrology. ISBN: 978-960-474-057-4 USCB (U.S. Census Bureau) 2014. TIGER/Line Roads Shapefile. United States Virgin Islands: https://www.census.gov/geo/maps-data/data/tiger-line.html U.S. Environmental Protection Agency. Guidance for Water Quality-based Decisions: The TMDL Process. April 1991. Prepared by Assessment and Watershed Protection Division U.S. Environmental Protection Agency. U.S. Environmental Protection Agency. Dissolved Oxygen TMDL for Mangrove Lagoon and Benner Bay, St. Thomas, U.S. Virgin Islands. September 15, 2003. Prepared for U.S. Environmental Protection Agency, Region 2. Prepared by U.S. Virgin Islands Department of Planning and Natural Resources. U.S. Environmental Protection Agency. Dissolved Oxygen TMDL for Vessup Bay and Red Hook Bay, St. Thomas, U.S. Virgin Islands. September 20, 2005. Prepared for U.S. Environmental Protection Agency, Region 2 and U.S. Virgin Islands Department of Planning and Natural Resources. U.S. Environmental Protection Agency. Dissolved Oxygen Total Maximum Daily Load, Salt River Bay, St. Croix, USVI. September 21, 2004. Prepared for U.S. Environmental Protection Agency, Region 2 and U.S. Virgin Islands Department of Planning and Natural Resources. U.S. Environmental Protection Agency. Dissolved Oxygen, Total Phosphorus, Turbidity, Fecal Coliform, and Enterococci Total Maximum Daily Loads, North Shore St. Croix Assessment Units, St. Croix, United States Virgin Islands. June 30, 2007. Prepared for: U.S. Environmental Protection Agency, Region 2 and U.S. Virgin Islands Department of Planning and Natural Resources. U.S. Environmental Protection Agency. Total Maximum Daily Load for St. Thomas Harbor, United States Virgin Islands. September 3, 2010. Prepared for U.S. Environmental Protection Agency, Region 2 and U.S. Virgin Islands Department of Planning and Natural Resources. USGS (United States Geological Survey). 2021. StreamStats: Streamflow Statistics and Spatial Analysis Tools for Water-Resources Applications. < https://www.usgs.gov/mission-areas/water- resources/science/streamstats-streamflow-statistics-and-spatial-analysis-tools>