w EPA United States Environmental Protection Agency velopmeril Watershed Classification Framework for the State of West Virginia L J. , *. / i ------- EPA/600/R-03/141 July 2004 Watershed Classification Framework for the State of West Virginia WV R-EMAP Final Report Naomi E. Detenbeck US EPA Mid-Continent Ecology Division National Health and Environmental Effects Research Laboratory Duluth, MN 55804 Leslie A. Jagger1 Stacey L. Stark1 OAO Corporation Duluth, MN 55804 Matthew A. Starry2 Computer Sciences Corporation Duluth, MN 55804 'EPA FAIR Contract, 2EPA FAIR II Contract CR82872001 Project Officer Naomi Detenbeck Mid-Continent Ecology Division National Health and Environmental Effects Research Laboratory Duluth, MN 55804 Mid-Continent Ecology Division National Health and Environmental Effects Research Laboratory Office of Research and Development U.S. Environmental Protection Agency Recyded/Recydable Printed with vegetable-based ink on paper that contains a minimum of 50% post-consumer fiber content processed chl ori ne free ------- Notice The information in this document has been funded wholly by the U.S. Environmental Protection Agency (US EPA), including support to U.S. Geological Survey EROS Data Center (IAG DW- 14938973) and West Virginia Department of Natural Resources (CR-82872001); and contract support to OAO Corporation under US EPA FAIR Contract 68-W5-0065, Delivery Order #24, and Computer Science (CSC) Corporation under US EPA FAIR II Contract 68-W01/W02-032, Task Order #024. This document has been prepared at the US EPA National Health and Environmental Effects Research Laboratory, Mid-Continent Ecology Division, in Duluth, Minnesota, with support from EPA FAIR II Contract 68-W01 /W02-032 to CSC. It has been subjected to review by the US EPA National Health and Environmental Effects Research Laboratory and approved for publication. Approval does not signify that the contents reflect the views of the Agency, nor does mention of trade names or commercial products constitute endorsement or recommendation for use. Preferred citation: Detenbeck, N.E., L.A. Jagger, S.L. Stark, and M.A. Starry. 2004. Watershed classification framework for the state of West Virgnia: WV R-EMAP Final Report. EPA/600/R-03/141. U.S. Environmental Protection Agency, Office of Research and Development, National Health and Environmental Effects Research Laboratory, Mid-Continent Ecology Division, Duluth, MN. ------- TABLE OF CONTENTS LIST OF ACRONYMS iv LIST OF FIGURES v LIST OF TABLES vii LIST OF APPENDICES vui INTRODUCTION I PURPOSE OF CLASSIFICATION AND INTENDED USE OF DATA 2 Purpose of watershed classification 2 Hydrology-based classification of watersheds 4 Relationship between HUCs and watersheds 5 Use of watershed classification as a sampling and assessment framework 6 METHODS 9 USGS watershed delineation 9 Delineation and coding of hydrologic units for the state of WV 9 Watershed characterization 14 Creation of sample frame and sample design for 12-digit HUCs 17 Response threshold derivation 20 CLASSIFICATION RESULTS 30 Hydrology thresholds 30 Distribution of land-use/land-cover variables across 12-digit subwatersheds in West Virginia 33 FUTURE UPDATES AND CLASSIFICATION REFINEMENTS 45 CONCLUSIONS 50 ACKNOWLEDGMENTS 51 REFERENCES 53 APPENDIX A: METADATA for MAPS AND DATABASES, INCLUDING DATA QUALITY INFORMATION A-l APPENDIX B: AMLS B-l APPENDIX C: DATABASE TABLES C-l 111 ------- LIST OF ACRONYMS DEM Digital elevation model DEP Department of Environmental Protection DNR Department of Natural Resources EDNA Elevation Dataset with National Applications EMAP Environmental Monitoring and Assessment Program FGDC Federal Geographic Data Committee GIS Geographic Information System MAIA Mid-Atlantic Integrated Assessment NHD National Hydrography Dataset NWBD National Watershed Boundary Database NWI National Wetlands Inventory PRISM Parameter-elevation Regressions on Independent Slopes Model R-EMAP Regional Environmental Monitoring and Assessment Program US EPA U.S. Environmental Protection Agency USGS U.S. Geological Survey i\ ------- LIST OF FIGURES Figure 1. Use of watershed classification 3 Figure 2. Definition of a subwatershed region 7 Figure 3. Generic example of Pfafstetter coding for an 8-digit HUC 8 Figure 4. USGS gaging stations and associated watershed boundaries 10 Figure 5. Map of 8-digit hydrologic cataloging units (HUCs), 10-digit watersheds, and 12- digit subwatershed units produced for West Virginia 11 Figure 6. Automated derivation of main channel slope 18 Figure 7. Major drainage basins associated with the state of West Virginia and surrounding states, as defined by 2- and 4-digit HUCs 19 Figure 8. Example of unequal weighting of interbasin HUCs 21 Figure 9. Map of ecoregions overlapping with West Virginia watersheds 22 Figure 10. Map of land-use covering West Virginia watersheds 27 Figure 11. Map of palustrine and lacustrine classes from National Wetlands Inventory for the state of West Virginia 29 Figure 12. Predicted 2-year flood normalized to watershed area as a function of watershed storage 31 Figure 13. Results of Classification and Regression Tree (CART) analysis 32 Figure 14. Map of WV 12-digit HUCs coded by fraction watershed storage 34 Figure 15. Map of WV 12-digit HUCs coded by estimated fraction impervious surface area .... 35 Figure 16. Map of WV 12-digit HUCs coded by estimated fraction surface mining area 36 Figure 17. Map of WV 12-digit HUCs coded by fraction agricultural area 37 Figure 18. Map of WV 12-digit HUCs coded by fraction forested area 38 Figure 19. Map of Omernik Level III ecoregions within West Virginia with land-use summaries. 40 Figure 20. Map of WV 12-digit subwatershed units by land-use and watershed storage cksses.. . 41 ------- Figure 21. Watershed cksses associated with 12-digit HUCs in sample population for Central Appakchkn Pkteau and Central Ridge and Valley ecoregions in West Virginia. ... 43 Figure 22. Watershed classes associated with 12-digit HUCs in sample population for Western Allegheny Plateau ecoregion in West Virginia 44 Figure 23. Effect of exponent term for storage on shape of relationship between watershed storage and normalized peak flow 48 Figure 24. Effect of peak flow breakpoint, exponent term, and prediction error on correct classification rate for high peak flow ckss or misckssification rate for low peak flow ckss 49 VI ------- LIST OF TABLES Table 1. Guidelines for hydrologic unit subdivision for National Watershed Boundary Database, according to interagency protocol (FGDC 2002) 13 Table 2. Geographic information system databases used to characterize WV watersheds 15 Table 3. Definition of sample frame for R-EMAP wadeable stream sampling in WV, 2000-2001 20 Table 4. Potential watershed characteristics related to peak flows 24 VII ------- LIST OF APPENDICES (CD-ROM, back pocket) APPENDIX A: INVENTORY OF MAPS AND DATABASES A-l APPENDIX B: AMLS B-l APPENDIX C: DATABASE TABLES C-l VIM ------- INTRODUCTION The Environmental Monitoring and Assessment Program (EMAP) is a research program to develop the tools necessary to monitor and assess the status and trends of national ecological resources. These tools include probability-based survey designs and indicators of biological condition. Regional EMAP (R-EMAP) is designed to evaluate how these tools can be applied at local and regional scales to meet the management needs of the states and regions. One of the emerging issues in monitoring programs for the states and tribes is the need to develop designs to meet multiple objectives outlined in different sections of the Clean Water Act. To meet the requirements of Section 305(b), probabilistic survey designs are needed to produce estimates of regional condition with a known level of confidence. In addition, monitoring programs must identify impaired waters and associated causes of impairment to meet the listing requirements of Section 303 (d) of the Clean Water Act (US EPA 2002). Ideally, monitoring programs will provide information to allow states to more efficiently plan the next round of monitoring, in order to identify as many impaired waters of the state as possible. One approach which can accomplish both of these needs is to incorporate risk-based categories into sampling designs, either as strata in a random-stratified design, or as categories to which different probability-weights are assigned. This approach allows managers to summarize not just information about regional condition, but also information about the risk of impairment associated with different classes of aquatic resources. A R-EMAP project was designed for the state of West Virginia to accomplish multiple objectives, including the evaluation of a risk-based random-stratified sampling process. An additional goal of the WV R-EMAP project was to refine the fish index of biotic integrity (IBI) developed through EPA's Mid-Atlantic Integrated Assessment (MALA) project to produce separate indices for different thermal classes of streams and specific to the state of West Virginia. The purpose of this report is to describe the methods involved in producing a preliminary watershed classification and incorporating the risk-based watershed classification into a monitoring design for ------- the state of West Vkginia. Validation and refinement of the watershed dassification scheme, based on analysis of pending monitoring data, will be addressed in a subsequent report. PURPOSE OF CLASSIFICATION AND INTENDED USE OF DATA Purpose of watershed classification Historically, classification systems have been used for inventory purposes, for stratifying landscapes prior to characterizing reference or regional condition, and for facilitating communication with natural resource managers and the public (Omernik and Gallant 1988; Heiskary and Wilson 1990; Herlihy et al. 2000; Pan et al. 2000; Waite et al. 2000). More recently, it has become necessary to develop classification systems that can explain differences in vulnerability of aquatic ecosystems to stressors, aid in regionalizing water quality criteria, and predict probability and causes of impairment of aquatic systems (US EPA 2002). Most sources of stream impairment contained in state 303(d) listings are related to nonpoint source pollution (US EPA 2001, 2003). To more efficiently deal with water quality management issues, an integrated approach to small watershed assessment, diagnosis, and restoration planning is needed (Figure 1). Current guidance from US EPA supports the development of a consolidated assessment and listing approach to enable joint application of Clean Water Act Sections 305(b) and 303(d) (US EPA 2003). Monitoring strategies developed for 305(b) regional assessments also must inform the 303(d) listing process. In the 2001-2003 WV Regional Environmental Monitoring and Assessment Project (R-EMAP), use of a new watershed classification technique specific for West Vkginia should allow the state to improve its ability to discern anthropogenic causes of impakment from non-anthropogenic sources of variability in reference condition (Cincotta 2000; Brazner et al. 2003) by producing a series of stressor-response relationships specific to different watershed classes. ------- threshold identification strategies: monitoring, modeling, restoration i I cost-benefit analyses watershed prioritization classification monitoring framework class condition and cumulative condition assessment screenm threshold confirmation [translators] Figure 1. Use of classification in sequence of monitoring, assessment of condition, and prioritization of watersheds for further monitoring, modeling, and restoration. In addition, use of the watershed classification scheme within a monitoring strategy will allow estimation of the condition of classes of watersheds and prediction of ecological risk for unmonitored systems. Gradients of land-use can be studied in combination with hydrology moderating factors to establish relationships of changes in runoff, and associated water quality and biological responses, with land-use activities (Verry 1986; Detenbeck et al. 2000). Since multiple non-point stressors potentially are affecting biological condition, it is appropriate to use watersheds as the sampling and/or assessment unit. In order to establish an appropriate stratified random sampling regime to diagnose causes of impairment, watersheds within the state of West Virginia must be properly classified. The classification developed and described in this document is based on identifying the levels of hydrology-moderating factors such as watershed storage and watershed development at which rapid degradation occurs. Watershed storage capacity is operationally defined as the fraction ------- of watershed area present in lake + wetland area, and can be used to predict the magnitude of floods (Jennings et al 1993). Peak flows are generally associated with a large fraction of nutrient and clean sediment yields from watersheds, so that watershed storage can also be used to predict attenuation of nonpoint-source loadings. Wetlands and lakes increase retention time within the watershed, allowing processes such as nutrient uptake, denitrification, and sedimentation in wetlands to improve water quality. Watershed storage can also be a good predictor of baseflow water quality (Johnston et al. 1990; Detenbeck et al. 2000; Detenbeck et al. 2003a,b). Hydrology-based classification of watersheds Watershed properties such as depressional storage volume and surface-water hydrogeomorphic types can be used to classify lakes or streams according to relative risk of impact. Our approach to predicting stream sensitivity to nonpoint source pollution is based on the nonlinear response of hydrologic regimes and associated loadings of non-point source pollutants to watershed properties (Richards 1990). Selected hydrologic thresholds are related to 1) variation in levels of watershed storage (either natural or anthropogenic) and 2) land-use activities affecting runoff and the hydrologic regime. In common usage, the term "threshold" refers to "the point at which a physiological or psychological effect begins to be produced" or "the point at which something starts." We define a hydrologic thresholds a breakpoint or inflection point in a nonlinear relationship between a watershed property and hydrologic response variable such as peak flows. For example, peak flows per unit area often increase exponentially once the storage capacity of a watershed has been exceeded. We have identified these thresholds based on a space-for-time substitution, or comparison across watersheds, as inadequate historical records are available for analysis of time series of individual watersheds. The United States Geological Survey (USGS) has defined a series of empirical nonlinear ------- equations relating watershed properties such as watershed area, channel slope, watershed storage, and land-use (percent forested, percent urbanization, or percent impervious surface area) to peak flows of given recurrence intervals (Q2, Q5,... Q100; Jennings et al. 1993). Peak flows increase exponentially as watershed storage decreases below a given threshold (Detenbeck et al. 2000). These thresholds form the basis for a watershed classification system. Relationship between HUCs and watersheds The first step in developing a classification strategy is to define the population of interest. Theoretically, an infinite number of watersheds could be defined, starting at any random point on a stream network (Detenbeck et al. 2003a). Benefits of defining a fixed set of watershed units include: a) the ability to tie assessment units to management plans and actions, b) ease in mapping classification units and classes, and c) the significance of fixed management units for stakeholders who form watershed management organizations to facilitate local environmental protection. The USGS is coordinating with other Federal agencies and the states to develop a National Watershed Boundary Database (NWDB), which will sequentially divide existing 8-digit Hydrologic Units (HUCs) into 10-digit HUCs, and 10-digit HUCs into 12-digit HUCs (Legleiter 2001). Use of a seamless national boundary database will facilitate communication among neighboring states, and planning and management of basins that cross interstate boundaries. In addition, creation of seamless nationwide GIS kyers will allow development of a series of elevation, hydrography, and derived coverages that are internally consistent with one another (Franken et al. 2001). We chose 12-digit HUCs as the base map for our classification framework. Watersheds associated with 12-digit HUCs span the size range of wadeable streams in West Virginia (400 - 40,000 ha), correspond to the size range of most management units in the state, and fit within a hierarchical scheme of larger management units (i.e., 10- and 8-digit HUCs). A total of 36 8-digit ------- HLJC units drain into the state of West Virginia; these 8-digit HUCs contain a total of 242 10-digit HUCs (USGS watershed units) and 1427 12-digit HUCs (USGS subwatershed units). Of these 12- digit HUCs, 883 fall predominandy within the borders of West Virginia. Not all HUCs are equivalent to watersheds. HUCs can be defined as either basin or interbasin. Noncoastal basin HUCs are generally equivalent with watersheds at the 12-digit scale, while interbasin HUCs divide the main channel of the next larger size class of HUCs into segments. Coastal HUCs can contain multiple parallel watersheds draining to a large water body (Great Lake or ocean), but these are not an issue in West Virginia. Interbasins may include a middle segment of a larger river that encompasses many smaller side tributaries. Interbasin HUCs can be aggregated with upstream units to define full watersheds (Figure 2). Throughout this report, derivation of watershed attributes for 12-digit HUCs is based on the boundary and characteristics of the full aggregated watershed regions. Use of Pfafstetter codes facilitates the discrimination between basin and interbasin HUCs, and identification of all upstream or downstream HUCs. Pfafstetter codes are not the same as HUC codes, but are produced by the USGS during the process of creating elevational derivatives such as watershed boundaries. Pfafstetter codes are hierarchical, with a 2-digit segment added at each new level of subdivision. The last digit of the Pfafstetter code increases consecutively from mouth to headwater sub-basins, with even numbers associated with tributaries (basin HUCs) and odd numbers associated with HUCs along the mainstem (interbasin HUCs) (Verdin and Verdin 1999; Figure 3). Use of watershed classification as a sampling and assessment framework Once a fixed set of watersheds is identified, delineated, characterized, and classified, this set can be used to define a sample population and survey design. At least two options exist for ------- Subwatershed Region HUC 050701010103 Region with National Hydrography Data and Shaded National Elevation Data Legend | | 050701010103 Region -/v NHD stream reaches Map Projection Alters Equal Area Conic Figure 2. Definition of a subwatershed region for 12-digit HUC 050701010103 through aggregation of 05 additional upstream HUCs. ------- Figure 3. Generic example of Pfafstetter coding for an 8-digit HUC (red outline). The base of the main channel is coded as "1" while the headwater portion of the main channel is coded as "9". Each of the four largest tributaries (basin HUCs) is assigned an even number between 2 and 8, in ascending sequence from the mouth towards the headwaters. Likewise, each of the segments draining to the mainstem between the tributaries (interbasin HUCs) is assigned an odd number between 3 and 7, in ascending sequence from the mouth to the headwaters. Coding of subdivisions within one of these 9 HUCs would proceed in a similar fashion, with the addition of a second digit to the code. ------- incorporating watershed classes into a survey design framework: 1) classes can be used as strata in a stratified random design, or 2) classes can be assigned unequal probability weights to influence the distribution of selected watershed units across classes (Detenbeck et al. 2003a). Both options allow a rare but significant watershed class (e.g., one with high associated probability of impact) to be sampled at a higher frequency than would be possible if a simple random survey design were implemented. Monitoring data can then be assessed for classes of watersheds as well as for an entire region, and a series of appropriate management strategies devised. METHODS USGS watershed delineation To define hydrologic thresholds, gaging stations with long-term flow records and derived peak flow statistics were identified from previous studies in the state of West Virginia (Frye and Runner 1970; Runner 1980; Figure 4). Watersheds associated with USGS gaging stations were delineated onscreen in Arclnfo using digital raster graphic (DRG) backdrops (1:24,000) and National Hydrography Dataset (NHD) coverages. Delineation and coding of hydrologic units for the state ofWV In contrast to manual delineations for USGS gaged watersheds, hydrologic units in West Virginia were delineated through a two or three stage process as part of the development of the Elevation Dataset with National Applications (EDNA, Franken et al. 2001; http://edcntsl2.cr.usgs.gov/ned-h/). The EDNA process produces HUC boundaries suitable for inclusion in the National Watershed Boundary Dataset (Figure 5). Delineation tasks were accomplished through an interagency agreement with the USGS EROS Data Center and collaborative efforts of US EPA Region III staff. Stage I of the EDNA process previously had been ------- USGS Gaging Station Watersheds Map Projection AJbar* Equal ATM Conic Gagng Slatnn Ports Gagng Station Walenhedi Level III Ecorrgions Name Wteslem Allegheny PbKau Central Appaladnam Cenlial Appalachun Ridgn and Va Figure 4. Location of USGS gaging stations and associated watershed boundaries used for derivation of hydrologic thresholds in West Virginia. 10 ------- WV8,10, &12-Digit HUCs with Contributing 8-Digit HUCs | J 8-Digit HUCs (EDNA) 8-Digit HUCs (contributing) 10-DigitHUCs 12-Digit HUCs Processed HUCs Data Type USGS1:250K EDNA Stage II EDNA Stage III Map Projection. Alters Equal Area Conic Figure 5. Map of 8-digit hydrologic catologing units (HUCs), 10-digit watersheds, and 12-digit subwatershed units produced for West Virginia. Stage III 12-digit HUCs are shown only for three of the 13 8-digit HUCs completed as of June 2003, i.e., those three 8-digit HUCs requiring additional processing beyond Stage II to correct flow directionality. 11 ------- accomplished through a memorandum of understanding between the USGS and US Weather Service. Stage I consists of an automated watershed delineation process using digital elevation models (DEMs) with GIS to produce a nationwide coverage of catchments at the scale of 1-2 mil2. Catchments produced at this stage reflect the predicted movement of surface water over surface topography; elements such as karst features are not included. Stage I produces several products, including catchment boundaries, synthetic streamlines (predicted location of streams and river courses), and elevational derivative raster coverages such as slope, aspect, flow direction, and flow accumulation. During Stage II, Stage I catchments are aggregated to the scale of 10- and 12-digit HUCs, based on a set of national guidelines (FGDC 2002; Table 1). At this point, GIS coverages are evaluated for internal consistency and potential errors are flagged. Indicators of potential errors either in digital elevation models or in NHD include divergence of synthetic streamlines from mapped NHD streams and inconsistencies between HUC boundaries and other sources of mapped watershed boundaries. During Stage III, a series of Arc View tools are used to correct small regions of the original digital elevation models near the flagged errors to represent the actual flow of surface water. For example, a highway overpass could be detected as a barrier to flow in an original DEM, when in fact, water can freely flow underneath the underpass. Creating a small "notch" in the original DEM at that neighborhood allows the stream course to be correctly predicted. After the hydrologically-corrected DEMs have been created, boundaries for 10- and 12-digit HUCs are regenerated, as are raster coverages for other hydrologic derivatives. 12 ------- Table 1. Guidelines for hydrologic unit subdivision for National Watershed Boundary Database, according to interagency protocol (FGDC 2002). Unit digits 8-digit 10-digit 12-digit Name Size range (acres) cataloging unit watershed 40j000 tQ 250>m subwatershed -, Q QQQ ^ QQQ Subdivisions/unit 5-15 5-15 For the West Virginia R-EMAP project, we used Stage II HUCs to develop the watershed classification framework. Currently, Stage III HUCs are only partially completed for the state, and were only used where Stage II HUCs required substantive correction, e.g., to correct direction of flow (Figure 5). The assignment of Pfafstetter codes during Stages I and II of the process allows ready identification of upstream and downstream 10- and 12-digit HUCs within 8-digit HUCs, and thus aggregation of interbasin HUCs to define watersheds using REGION processes in Arclnfo. Definition of Arclnfo regions is desirable considering the nested nature of watersheds. Regions are model areal geographic features described by one or more polygons. Regions are efficient ways to model nested or overlapping data because of the concept of shared geometry. For example, one polygon feature can belong to multiple regions without the topology being duplicated. Two files store the region-arc relationship and the region-polygon relationship "behind the scenes." Regions are stored as subclasses within a polygon coverage and each region subclass has its own separate attribute table. The bulk of the work was completed using individual Arc Macro Language (AML) scripts developed for each component of the characterization process. The scripts were used to create and combine region subclasses into an integrated coverage for further geoprocessing using the 13 ------- REGIONQUERY command in Arclnfo (ESRI, Redknds CA). The REGIONQUERY command mimics traditional oveday commands on region subclasses in an integrated coverage. So, for example, if a single coverage had a 12-digit watershed subclass and a wetland subclass, the REGIONQUERY command could be used to create a new subclass computed from the geometric intersection of the two subclasses (watersheds and wetlands). The region tables keep track of which wetland polygons belong to each watershed polygon while only physically representing the polygons once. One table could then be created that summarized wetland type by area for each watershed region. Due to the size of input datasets and software/hardware limitations, 8-digit HUCs containing 12-digit watersheds were processed individually. This presented an issue when a 12-digit watershed required additional upstream 8-digit HUC(s) to be considered a true watershed. A method was needed to identify diese "interbasin" watersheds and combine the output for their required counterparts. Each interbasin watershed was given a special identifier code in the watershed region attribute table. A table was made listing all possible interbasin codes and the 12- digit watersheds needed to complete the full watershed. This table was imported into a Microsoft Access database. The individual output tables for 12-digit watersheds within an 8-digit HUC were combined into one table for each characterization process and imported into the same database. A series of queries were performed to combine the individual output values to determine the total output values for interbasin watersheds for each characterization process. Watershed characterisation Watershed regions associated with all 12-digit HUCs in the state of West Virginia and long- term gaging stations were characterized for attributes expected to be related to generation of peak flows (Fne and Runner 1970; Runner 1980): watershed area, main channel length, main channel 14 ------- slope, watershed storage, mean January minimum temperature, snowpack, and percent forest cover. In addition, watersheds were characterized for other major land-use and land-cover characteristics related to generation of runoff and nonpoint source pollution: percent mined land, percent impervious surface area, and percent agriculture. GIS databases used to calculate watershed characteristics are described in Table 2. Table 2. Geographic information system databases used to characterize West Virginia watersheds. EDNA = Elevation Dataset with National Applications, NHD = National Hydrography Database, NWI = National Wetlands Inventory, OWI = Ohio Wetlands Inventory, NYWI = New York Wetlands Inventory, PRISM = Parameter-elevation Regressions on Independent Slopes Model. Database EDNA NHD NLCD- mining modified NWI OWI PRISM Layers and attributes 10-digit watershed boundaries 12-digit watershed boundaries Pfafstetter codes elevation grids stream reaches stream level land-use, with updates for surface mining palustrine wetland polygons lacustrine polygons wetland and open water types by grid cell monthly snowfall annual precipitation minimum January temperature Derived variables sample units watershed boundaries watershed areas upstream connections main channel slope main channel slope upstream connections percent land-cover percent wetland + lake area percent wetland + lake area snowfall mean annual precipitation mean minimum January temperature Last update Scale Source 2003 1:24,000 USGS EROS Data ,.. ... Center (30 m gnd) 1999 1:100,000 USGS 2001 1:24,000 NLCD - EPA (30 m grid) surface mining grid - Tennessee Valley Authority Thru: US EPA Region 3, Wheeling, WV 1981 to 1:24,000 US FWS, WV DEP present composite coverage based on 1:24,000 Ohio Department of 1985-1987 Natural Resources TM scenes 1998, based 1:250,000 Climate Data Source, on 1961-1990 Corvallis, OR climate data, 15 ------- The PRISM (Parameter-elevation Regressions on Independent Slopes Model; Climate Source, Corvallis, OR) GIS coverage was used to generate climatic summaries across watersheds. PRISM is a continuous grid (raster coverage) of climatic variables generated from weather station data and regional nonlinear regressions that take into account regional precipitation averages as well as orographic effects associated with altitude, and proximity/influence of large surface water bodies (Daly et all 994). We used National Wetlands Inventory (NWI) coverages, supplemented by the Ohio Wetlands Inventory coverages to define watershed storage. Total areal coverage of palustrine and lacustrine polygons was added and divided by watershed area. For the Ohio Wetlands Inventory coverage, polygons coded as 34 to 38 (wetland and open water cksses) were included. For the New York Wetlands Inventory, cover classes 1-4 and 9 were included (all classified or unclassified wetland polygons). Fraction impervious surface area was derived as a function of knd-use intensity as: fimperv = (0.55*flwintrs)+(0.90*£hintres)+fcomind, where fimperv = fraction impervious area in watershed, flwintrs = fraction low intensity residential area in watershed, fhintres = fraction high intensity residential area in watershed, and fcomind = fraction commercial, industrial, and transportation knd-use area in watershed. based on median of estimated constructed materials associated with different land-use cksses in NLCD (http://www.epa.gov/mrlc/dennitions.html). These weightings are higher than those 16 ------- reported for northern Virginia (NVPDC 1980), because the latter included only effective impervious surface area (unconnected rooftops were not included), but are consistent with classification guidelines for NLCD. Main channel slope was defined first by selecting main channel reaches from the NHD coverage based on the minimum value for the reach LEVEL attribute within a watershed. Main channel reaches were buffered by 30 meters, and this buffer was used to clip a subset of elevation grids from the DEM coverages. Frequency analyses were used to estimate the 10th and 85th percentile of elevation values from the stream elevation buffer files; combined with main channel length, these could be used to calculate main channel slope (Figure 6). Creation of sample frame and sample design for 12-digit HUCs Although we developed the classification framework for the entire state of West Virginia, in practice, we restricted membership in the sample frame for the 2001-03 WV R-EMAP project. Sample frames were defined based on drainage basin, ecoregion, and 12-digit HUC watershed region size (Table 3). The Potomac River basin was excluded from study because of species differences related to biogeography, and because the fish IBIs developed for adjacent regions in Maryland should be adequate for the state of WV to apply in that basin (Figure 7). In practice, we also excluded watersheds with ill-defined drainage networks related to karst topography (e.g., no single outlet, underground streams). Selection of 12-digit HUCs from within the sample frame each year was based on a random-stratified design. Strata were based on watershed classes defined by hydrologic thresholds and knd-use intensity. After defining 12-digit HUCs as single independent watersheds (basin HUCs) or groups of hydrologically adjacent (dependent) interbasin units, 17 ------- HUC 060701010103 Region with NHD and Shaded NED / c "Ajf-WI. 1 0 Ot5 01 N A MUC 050701010103 Region wtm NMD nd ShxtedNED 0.8 0.6 0.4 0.2 Derivation of main channel slope 85% 10% 472 - 436 = 36 m / change in elevation 36 m elev/ 29.928 km main channel length 0 420 430 440 450 460 470 480 490 500 510 Elevation (meters above sea level) Figure 6. Automated derivation of main channel slope. The main channel is identified based on finding stream reaches with the minimum "LEVEL" attribute in NHD (top left), then is buffered and intersected with a digital elevation model (top right). A frequency analysis of output elevation grid cells yields estimates of the 10th and 85lh percentiles. Combined with main channel length, these can be used to calculate main channel slope. ------- Major Drainage Basins Map Projection: Alters Equal Area Conic | | 2-Digit HUCs | | 4-Digit HUCs HUC 02 02 05 HUC 04 0207 0208 0501 0502 0503 0504 0505 0507 0509 REGION NAME Mid Atlantic Region Ohio Region SUBREGION NAME Potomac Lower Chesapeake Allegheny Monongahela Upper Ohio Muskingum Kanawha Big Sandy-Guyandotte Middle Ohio Figure 7. Major drainage basins associated with the state of West Virginia and surrounding states, as defined by 2- and 4-digit HUCs. 19 ------- unequal probability-weighting was used to limit the chance of selecting subwatersheds directly upstream or downstream of one another (Figure 8). Table 3. Definition of sample frame for R-EMAP wadeable stream sampling in West Virginia, 2000-2001. Year Ecoregion(s) Drainage basin Watershed size (ha) 2000 Central Appalachian Pkteau inclusive 400-40000 2000 Central Ridge and Valley all excluding Potomac River 400-40000 Ecoregion basin 2001 Western Allegheny Pkteau inclusive 400 - 40000 Response threshold derivation We carried out analyses and applied sampling designs separately by ecoregion for studies in 2000 and 2001 (Table 3, Figure 9). To define hydrologic thresholds related to natural watershed attributes, we identified gaging stations with long-term flow records and derived peak flow statistics from previous studies in the state of West Virginia (Frye and Runner 1970; Runner 1980; Figure 4). We delineated watersheds for each gaging station and characterized these watersheds based on attributes recorded in Frye and Runner (1970) and Runner (1980). We updated watershed storage estimates from Frye and Runner's values using GIS analysis of National Wetlands Inventory coverages and calculation of fraction watershed area covered by palustrine and kcustrine systems. In earlier analyses conducted before GIS was widely avaikble, Frye and Runner calcukted watershed storage by overkying a point grid on USGS topographic quadrangles and counting intersections with wetlands, lakes, and ponds. With a rare aquatic resource such as wetknds occurring in West Virginia or other regions with well-developed drainage systems, the accuracy of such methods is limited by the temporal accuracy of maps, the density of grid points used, and number of grid 20 ------- Kilometers 25 50 100 A Random-Stratified Selection Process Weighting Procedure Map Projection Alters Equal Area Conic Legend | | Independent 12-digitHUCs | | Dependent 12-digit HUCs i Central Appalachians Inclusion Probabilities | 0.013 0.018 0026 0.053 Figure 8. Example of unequal weighting of interbasin HUCs for one watershed class in the Central Appalachian Plateau ecoregion of West Virginia. 21 ------- Map Pro)*cton Albert Equal Aiea Cone US EPA Level III Ecoregions Level III Ecoregions Name Wo3»m Mogfwny Fbt.au C«n0.il AppatichKuij Cones! AppaLxrhoo Rrioo, and Vifay^ Figure 9. Map of ecoregions overlapping with West Virginia watersheds. 22 ------- points counted per unit area. For West Virginia, percent coverage had only been recorded to a single significant unit (as 0 or 1). We applied two different methods to identify hydrologic thresholds, i.e., regions in plots of area-normalized peak discharge (e.g., Q2/watershed area) showing nonlinear discrete shifts along a gradient of watershed attributes. We obtained values for 2-year peak flows from Frye and Runner (1970) and Runner (1980). We chose a recurrence interval of two years as an endpoint because this is the frequency of flooding associated with development of channel morphology and sediment delivery (Rosgen 1996). Our first method, applied in 2000, involved screening potential parameters of interest identified by Frye and Runner (1970; see Table 4) using Linear regression analysis of log- transformed Q2/area, with observations weighted by period of record to account for differences in level of uncertainty (Tasker and Stedinger 1989). We chose the best number and combination of predictive variables using Mallow's Cp statistic (SAS 1990). We identified thresholds through visual graphical analysis, by plotting (nontransformed) Q2/area as a function of single variables or combinations of variables and observing where nonlinearities in relationships were observed. In 2001, we applied a second technique, based on classification and regression tree (CART) analysis of potential variables in SYSTAT (Wilkinson 1999). Nonlinear regression analysis followed by visual analysis of graphic bi-plots was less successful in 2001 because of collinearities in independent variables. CART is a nonparametric technique that sequentially bisects observations into subpopulations based on a single response variable, and identifies dependent variables associated with breaks in the distribution. In this way, CART is able to identify nonlinearities in relationships, as well as interactions among predictor variables, based on a minimum set of assumptions (Wilkinson 1999). 23 ------- Table 4. Potential watershed characteristics related to peak flows, Frye and Runner (1970). Code A S St F L E S, Units mil2 feet -mil"1 Si T, Definition watershed area channel slope storage (area in lakes and ponds +1) forest area main channel length mean basin elevation snowfall index average annual precipitation - 20 precipitation intensity/24 hours, 2-year inches recurrence interval soil infiltration index average minimum January temperature deg F miles feet above sea level inches inches We identified potential thresholds of response related to land-use change from the literature because insufficient data were available from gaged sites with significant periods of record. We based potential land-use thresholds of response on previous analyses of data from the Maryland Biological Stream Survey (MBSS, Boward et al. 1999). For the MBSS data set, species loss in fish communities had been found at levels of impervious surface area above 2% of watershed area. Nonlinear shifts in water quality were associated with levels of agricultural land-use above 25%. Other studies have shown a change in peak snowmelt associated with loss of mature forest cover greater than 40-60% (Verry 1986). This approach is similar to trying to bracket expected response diresholds in bioassays for suspected toxicants. Less information is available from the literature to detect thresholds of impact rekted to mining activity in the Appalachian region. Scott (1984) identified percent disturbed area as the best 24 ------- predictor of changes in stream hydrology in a subset of mined watersheds of West Virginia, presumably associated with changes in water quality (Scott 1984). However, Scott did not have long hydrologic time series for these watersheds, but instead based his conclusions on modeled results; his analysis showed an apparent threshold at 30% disturbed area. We chose to use a conservative estimate of 10% for the mining threshold, similar to that shown as a threshold for percent impervious surface area in most literature reviews (Schueler 1994), with the knowledge that the magnitude of mining thresholds would have to be refined in the future. We defined land-use coverage for watersheds associated with 12-digit HUCs by intersecting watershed boundaries with a mining-modified NLCD coverage (Table 2, Figure 10). We defined percent watershed storage by intersecting watershed boundaries with wetland inventory coverages (Figure 11). 25 ------- Mining-Modified NLCD Land-Use/Land-Cover Mining-Modified NLCD Grouped Land-Use Codes Mining Agriculture/Orchards Forest Grassland Urban/Residential Open Water/Wetlands Barren Map Projection: Albers Equal Area Conic Figure 10. Map of land-use covering West Virginia watersheds, based on National Landcover Characterization Database, updated for mining category. 27 ------- 0 25 50 Map Projection: Abets Equal Area Conic Lakes and Wetlands in National Wetland Inventory (NWI) Legend 12-digit HUCs NWI SYSTEM Lacustrine Palustrine Figure 11. Map of palustrine and lacustrine classes from National Wetlands Inventory for the state of West Virginia. 29 ------- CLASSIFICATION RESULTS Hydrology thresholds In year 2000, the nonlinear equation we derived to describe 2-year peak flows (Q2) for the Central Appalachian Plateau and Central Ridge and Valley ecoregions, was: Q2 = Aa Pp T Snsn St5' (p < 0.05) where Q2 2-year peak flow A watershed area P = average annual precipitation T = average minimum January temperature Sn = snowfall St = watershed storage. Percent forest was not included in the equation for 2-year peak flow, probably because of a lack of significant variation in percent forest cover in these two ecoregions. When gaging station data were combined for the full state, percent forest cover was included as a significant predictor variable in peak flow equations for 5- and 10-year recurrence intervals, but not for 2-year recurrence intervals. Visual graphical analysis showed that percent forest alone showed no threshold effect on peak flows normalized to watershed area, but the product of percent forest cover, snowfall, and minimum January temperature did show a noisy threshold with a lot of scatter. The sharpest threshold for Q2/A identified based on visual graphical analysis was for percent watershed storage (Figure 12), with a threshold identified at ~0.5% storage (fraction storage = 0.005). We identified variables and associated threshold levels for Q2/A for the Western Allegheny Plateau in 2001 through CART analysis. At each step, the CART procedure sequentially divided observations into two categories based on the magnitude or category of the predictor variable (in 30 ------- 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Fraction storage (lake + wetland/wshd) Figure 12. Predicted 2-year flood normalized to watershed area (cfs/mi!2) as a function of watershed storage (lake + wetland area/watershed area) for USGS gaging station watersheds in the Central Appalachian Plateau and Central Ridge and Valley ecoregions of West Virginia. this case fraction storage or main channel length) that best separated the response variable, Q2/watershed area, into low and high magnitude subpopulations. Each box in the output is basically a frequency plot, with each observation represented by a colored dot, and each color representing the final subclass assignment. Separation of the dataset based on combination of two variables accomplished a 67% reduction in total variance (Figure 13). Flashy hydrologic regimes (significantly larger Q2/A) for the Western Allegheny Plateau ecoregion occurred in low storage watersheds (< 0.03% watershed storage) or relatively short drainage basins (main channel length < 13.5 km). 31 ------- Q2/watershed area I %storage<0.031 Main channel length (km)<13.5 Figure 13. Results of Classification and Regression Tree (CART) analysis of 2-year flood normalized to watershed area (cfs/mil2) as a function of watershed attributes for USGS gaging station watersheds in Western Allegheny Plateau ecoregion of West Virginia (percent reduction in error = 67 %). At each step, the CART procedure sequentially divided observations into two categories based on the magnitude or category of a predictor variable (in this case fraction storage or main channel length) that best separated the response variable, Q2/watershed area, into low and high magnitude subpopulations. Each box represents a frequency plot, with each observation represented by a colored dot, and each color representing the final subclass assignment. ------- Distribution of land-use/ land-cover variables across 12-digit subwatersheds in West Virginia Overall, the level of watershed storage in West Virginia is very low. For all watersheds associated with 12-digit HUCs in the 8-digit HUCs characterized (those within or overlapping with WV state borders), the top 20th percentile of HUCs had a range of 0.8 to 18.6% watershed storage (Figure 14). Most 12-digit HUCs in the top 20th percentile were actually located near the border or even outside of WV state boundaries, in the headwaters of the Upper and Middle Ohio River basins. Nonetheless storage thresholds associated with a marked change in hydrologic regime were low enough so that both high and low storage classes could be identified within each set of ecoregions (Figure 14). The range of percent impervious surface area in watersheds associated with the 12-digit HUCs characterized is also relatively low (0-18.4 %), with the top 20th percentile covering a range of only 1.2 - 18.4 %. The spatial distribution of the top 20th percentile of percent impervious surface area is similar to that of the top 20th percentile of watershed storage, paralleling development along road or river corridors, and probably the creation of reservoirs (Figure 15). As for urban development, the highest concentration of agriculture within West Virginia drainage basins (top 20th percentile = 48-83 %) occurs outside of the state boundaries in the Upper Ohio or Potomac drainages (Figure 16). In contrast, the highest concentration of mining activity (top 20th percentile = 1.4-16.4% watershed area) occurs within the Central Appalachian Plateau, in the southwestern portion of the state (Figure 17). Forest cover across the state is relatively high and constant; most instate HUCs are in the top 40th percentile of HUCs characterized (79 - 100 % forest cover; Figure 18.) 33 ------- Fraction Watershed Storage By 12-DigitHUC 12-Digrt HUCs FractionStorage 0.0000-0.0010 00011 -00023 00024 - 0.0042 HH 0 0043 - 0 0082 00083-0.1860 klip Protection Albett tqujl Aiea Come Figure 14. Map of 12-digit HUCs within 8-digit HUCs overlapping with West Virginia, coded by fraction watershed storage ([lake + wetland area/watershed area] x 100). Classes mapped are not equal interval classes, but represent population quintiles. ------- Fraction Impervious Surface By12-DigitHUC 12-Digit HUCs Fraction Impervious Surface 0.0000 - 0.0008 0.0009 - 0.0028 0.0029 - 0.0053 0.0054-0.0121 0.0122-0.1837 Map Projection. Abers Equal Area Conic Figure 15. Map of 12-digit HUCs within 8-digit HUCs overlapping with West Virginia, coded by estimated fraction impervious surface area in associated watersheds. Classes mapped are not equal interval classes, but represent population quintiles. 35 ------- Fraction Mining By12-DigitHUC 12-Digit HUCs Fraction Mining 0.0000 00001-0.0012 00013-0.0064 IB} 0.0065-0.0138 Ij^B 0.0139-0.1644 Map Pto)*ctKm *Jb«f* Equal Aiea Cone Figure 16. Map of 12-digit HUCs within 8-digit HUCs overlapping with West Virginia, coded by estimated fraction surface mining area in associated watersheds. Classes mapped are not equal interval classes, but represent population quintiles. 36 ------- Fraction Agriculture By 12-DigitHUC 12-Digit HUCs Fraction Agriculture 0.0002 - 0.0822 0.0823-0.1910 0.1911 -0.3140 . 0.3141-0.4814 ^^1 0.4815-0.8315 Map Projection: AJbers Equal Area Conic Figure 17. Map of 12-digit HUCs within those 8-digit HUCs overlapping with West Virginia, coded by fraction agricultural area in associated watersheds. Classes mapped are not equal interval classes but represent population quintiles. 37 ------- Fraction Forest By 12-DigitHUC 12-Digrt HUCs Fraction Forest 0.1400 - 0.4788 0 4789 - 0.6605 0 6606 - 0 7870 IB 07871 -08937 ^H 0 8938 - 0.9989 M«p Pioftctlon Albef* Equal Ai«« Come Figure 18. Map of 12-digit HUCs within those 8-digit HUCs overlapping with West Virginia, coded by fraction forested area in associated watersheds. Classes mapped are not equal interval classes but represent population quintiles. 38 ------- Because of the uneven distribution of watershed storage, watershed morphology, and land- use across the state of West Virginia, the distribution of 12-digit HUCs is not even across watershed classes; in fact some combinations of watershed classes are not even represented (Figures 19-20). For example, in the Western Allegheny Plateau ecoregion, only one 12-digit HUC watershed occurred in the high urbanization/low storage class, and no 12-digit HUC watersheds occurred with combinations of moderate or high agriculture and low storage. Thus, not all potential watershed classes can be compared within each ecoregion (Figures 21-22). 39 ------- Mdp Pto|*t.t)on Alb*ft Equal Ai»« Cone Aggregated Land-Use by Ecoregion [ | HUC Bomdwy Agriculture/Orchardi BH Forest Grasslands Urban/fteskfartfal H Open Water/Wetlands Barren Level III Ecoregions Name Western Allegheny Plateau Central Appalachians Central App. Ridges & Valeys Figure 19. Map of Omemik Level III ecoregions within state of West Virginia, and associated land-use summaries at ecoregion level for region within WV watersheds. 41) ------- Western Allegheny Plateau Central Appalachian Ridges & Valleys 12-DigitHUC Watershed Classes Watershed Classes Low Impact/Low Storage j Low Impact/High Storage High Urbanization/Low Storage High Urbanization/High Storage Moderate Agriculture/Low Storage Moderate Agriculture/High Storage High Agriculture/Low Storage High Agriculture/High Storage Ecoregions Ecoregions Level Level Map Projection: Albers Equal Area Conic Figure 20. Map of 12-digit subwatershed units within West Virginia's 8-digit HUCs, classified by land-use and watershed storage classes. Storage classes are based on empirical hydrologic thresholds for watershed storage (CA, CARV ecoregions) or combinations of storage and main channel length (WAP ecoregion). 41 ------- Map Projection: Abere Equal Area Conic Random-Stratified Selection of HUCs for Sampling Legend | | HUCs selected for sampling Thresholds Reference, low storage HI Reference, high storage High mining, low storage High mining, high storage BH High agriculture, low storage High agriculture, high storage High impervious, low storage i | High impervious, high storage Level III Ecoregions Name | | Central Appalachians Figure 21. Watershed classes associated with 12-digit HUCs in sample population for Central Appalachian Plateau and Central Ridge and Valley ecoregions in West Virginia. Sample population included 12-digit subwatersheds with associated watersheds in size range of wadeable streams and outside of Potomac River basin. 43 ------- Watershed Classes in Western Allegheny Plateau Map Projection. A!b« r_ Equal Area Conic Legend Sample Ponts Watershed Classes Low Impact/Low Storage Low Impact/High Storage | High Urbanization/Low Storage High Urbanization/High Storage Moderate Agriculture/Low Storage Moderate Agncuttvjre/High Storage | High Agriculture/Low Storage High Agriculture/High Storage Level III Ecoreglons Name 1 Western Allegheny Plateau Figure 22. Watershed classes associated with 12-digit HUCs in sample population for Western Allegheny Plateau ecoregion in West Virginia. Sample population included 12-digit HUCs with associated watersheds in size range of wadeable streams. 44 ------- FUTURE UPDATES AND CLASSIFICATION REFINEMENTS Potential future improvements to the watershed classification system for West Virginia include the following: Examination and documentation of differences in reference condition among low and high storage reference ("low impact") classes based on WV R-EMAP data, EMAP data from the MidAtlantic Integrated Assessment (MAIA), and macroinvertebrate and water quality data collected as part of the WV Department of Environmental Protection (DEP) rotating basin monitoring program Examination and documentation of differences in condition between reference and "impacted" watershed classes within storage categories based on WV R-EMAP data, EMAP MAIA data, and macroinvertebrate and water quality data collected as part of the WV DEP rotating basin monitoring program Examination and adjustment of thresholds of impact associated with watershed storage and knd-use intensity based on the analyses described above. Potential refinement of watershed impact classes based on empirical analysis of existing data with the addition of other watershed attributes such as development within stream and river corridor buffers, road density, and density of stream or river road crossings. Land-use impacts in a mountainous region might be better reflected by estimation of valley development, as compared to land-use intensity on a whole-watershed basis. Improved estimation of percent impervious surface area in WV watersheds. The method used for estimation in this project was indirect. However, the Chesapeake Bay program is currently mapping impervious surface area directly based on satellite imagery from Landsat 7, IKONOS, and MODIS platforms (http://chesapeake.towson.edu/data/), and 45 ------- additional studies are underway to calibrate percent impervious surface area estimation equations based on actual measurements from aerial photographs (David Jennings, US EPA, Reston, VA, personal communication). Improved estimation of percent mined area (surficial disturbance). Current estimates were based on interpretation of TM + SPOT satellite imagery by Tennessee Valley Authority under contract to Canaan Valley Institute for permitted mining areas (Hope Childers, US EPA Region III). The state of West Vkginia is currently constructing coverages of mining activity that might eventually allow better definition of area! extent and different sources of mining impacts (e.g., valley fill vs mountain top mining, http://www.nrac.wvu.edu/hollowfill/, http://gis.dep.state.wv.us/data/omr.html); however, current coverages are based on permit boundaries which do not always reflect mining activity. Slight adjustments in boundaries and characterizations from Stage II to Stage III 10- and 12-digit HUCs as Stage III EDNA process is completed for the state of West Vkginia. There are several ways in which derivation of thresholds could be improved, depending on the objectives of a monitoring project. CART analysis objectively maximizes the separation of two subpopulations of a dependent variable (e.g., Q2/A) as a function of the best predictive variable, e.g., fraction watershed storage. Thresholds derived from graphical analysis could be verified through application of a statistical procedure such as CART (Wilkinson 1999). Reanalysis of USGS data for the Central Appalachian ecoregions using CART yielded an apparent threshold of fraction storage = 0.001, or 0.1%, very close to the threshold of 0.5% derived through visual examination). 46 ------- Alternatively, thresholds could be chosen by optimizing the chance of correctly identifying watersheds with high normalized peak flows, i.e., those with greater potential for bank erosion and transport of nonpoint source pollutants, while minimizing the chance of misidentifying watersheds expected to have low normalized peak flows. If watershed storage is used as the indicator of normalized peak flows, then the probability of the two types of misclassification error will depend on the strength of the relationship between normalized peak flows and watershed storage, the error associated with the relationship between the two variables, and the underlying distribution of watersheds across the range of watershed storage. The absolute value of the exponent term in the equation relating storage and peak flows is a measure of the "strength" of the relationship. An exponent of zero would indicate that there is no relationship between storage and normalized peak flows, while a highly negative value indicates a rapid rate of increase in peak flows as storage is reduced (Figure 23). Several examples are shown below to illustrate the combined effect of the magnitude of the exponent (-2 versus -0.01), the magnitude of a desired breakpoint in Q/A relative to its underlying distribution (median versus 95th percentile), and the effect of differing levels of error of prediction (standard error of estimate = 0 to 100%; Figures 24a-h). In these simulations, an even distribution of watersheds across storage cksses was used. In all cases there is a tradeoff between maximizing the proportion of correct predictions for high normalized peak flows (Figures 24a,c,e,g) and minimizing the proportion of incorrect predictions for the low normalized peak flow class (Figures 24b,d,f,h). As we reduce the chosen threshold for percent watershed storage, our ability to predict watersheds with high normalized peak flows improves, but the misclassification error rate for the low peak flow ckss also increases (e.g., compare Figure 24a with Figure 24b). Our chances for 47 ------- Effect of exponent term for storage STA-5 STA-2 STA-1 STA-0.5 STA-0.2 STA-0.1 STA-0.01 STA0 20 40 60 80 100 120 % watershed storage Figure 23. Effect of magnitude of exponent term for storage on shape of relationship between watershed storage and normalized peak flow (Q/area = ST). making correct high flow predictions and minimizing errors in predicting low flow watersheds arc generally better the stronger the relationship is between storage and peak flows (e.g., compare Figure 24a with Figure 24c and Figure 24b with Figure 24d). The greater the error in the relationship between storage and normalized peak flows (higher SEE), the worse our predictions will be. (Compare different color lines within Figure 24a or 24b.) Flowever, prediction error appears to have little influence if we wish to distinguish only those watersheds with the most extreme peak flows (> 95lh percentile of Q/A; Figures 24e and f). In this example, our ability to avoid making misclassificanon errors for low peak-flow watersheds is relatively small if our breakpoint for Q/A is high (Figures 24f and h). If at the same time, the strength of the relationship is weak (e.g., exponent = -O.U1), however, we also have little chance of correctly identifying the high peak-flow watersheds (l-igurc 24g). 48 ------- a) c) e) g) Fraction LO storage dass with Q/A > median . s s e s - s Eft I c< 1 in 4 « | 1 i» 1" | 0.4 S 0.2 I o S. Effect of standard error of estimate on correct prediction of Q/A > threshold in LO storage class, exponent = -2 K^V^ - SEE = 0% - SEE = 1% - SEE=10% SEE=50% - SEE=7S% - SEE=100% > 20 40 SO 80 100 120 Chosen threshold, % watershed storage ect of standard error of estimate on >rrect prediction of Q/A > threshold LO storage class, exponent = -0.01 )hf^~~^ ^ SEE=60% - SEE=7S% - SEE=100% 0 20 40 60 80 100 120 Chosen threshold, % watershed storage | Eff £ CC < in a « 1 i 1« r §04 I" o fo | Eff 1 « » in i« 1 ' §0.8 D |« £0.4 2 0.2 1 ect of standard error of estimate on irrect prediction of Q/A > threshold LO storage classes, exponent = -2 V - SEE = 0% - SEE = 1% - SEE=10% SEE=SO% - SEE=75% - SEE=100% 9 20 40 60 80 100 120 Chosen mreshold.% watershed storage ect of standard error of estimate on rrect prediction of Q/A > threshold LO storage class, exponent = -0.01 SEE=SO% - SEE=76% - SEE=100% 1 ° £ 0 20 40 60 SO 100 120 Chosen threshold, % watershed storage b) d) f) Effect of standard error of estimate on c incorrect inclusion of Q/A > threshold I in HI storage class, exponent = -2 § 0.7 f 0.6 S 0.6 I** |0.3 V 0.2 V - SEE = 0% - SEE = 1% - SEE=10% SEE=60% - SEE=76% - SEE=100% £ 0 20 40 60 60 100 120 Chosen threshold, % watershed storage Eff c in< TO . 1 ln f 0.8 5 °-7 |0.6 S 0.6 CD « °A i 0.2 ts o C ect of standard error of estimate on correct inclusion of Q/A > threshold HI storage class, exponent = -0.01 __^A SEE=60% - SEE=75% - SEE=100% 9 20 40 60 80 100 120 Chosen threshold, % watershed storage S Eff 8 in( 1 0.07 £ 0.06 1 0.05 § 0.04 | 0.03 £ 0.02 * 0.01 ect of standard error of estimate on correct inclusion of Q/A > threshold n HI storage class, exponent = -2 I - SEE = 0% - SEE = 1% - SEE=10% SEE=SO% - SEE=7S% - SEE=100% 2 0 20 40 60 80 100 120 Chosen threshold, % watershed storage . Efl g in ? in < 0.07 1 0.06 S 0.04 f 0.03 a 0.02 S 0.01 I ° ect of standard error of estimate on correct inclusion of Q/A > threshold HI storage class, exponent = -0.01 /V\A SEE=50% - SEE=76% - SEE=100% 0 20 40 60 80 100 120 Chosen threshold, % watershed storage Figure 24. Effect of peak flow breakpoint (Q/A > median [a-d] versus Q/A >95th percentile [e-h]), magnitude of exponent term in relationship, Q/area = STa (a = -2 [a,b,e,fj versus a = -0.01[c,d,g,h]) and prediction error on correct classification rate for high peak flow class [a,c,e,g] or misclassification rate for low peak flow class [b,d,f,h]. 49 ------- CONCLUSIONS In some regions such as West Vkginia, the watershed classification framework for monitoring will be limited initially by available regional data and information on relationships between land-use and hydrologic or biological response. However, the process is iterative and classification schemes can be validated and improved upon based on initial monitoring results. In previous approaches to watershed classification, we have developed thresholds based on land-use attributes calculated as a fraction of watershed area. Land-use variables described as a fraction of watershed area might be less useful in montane regions where development is concentrated in valleys; in these cases it might be more appropriate to express land-use as a fraction of stream buffer zone area. Watershed storage was one to two orders of magnitude lower in West Vkginia than in Great Lake regions for which we have derived hydrologic thresholds for storage. However, the hydrologic threshold for watershed storage was also an order of magnitude lower for WV as compared to the upper Midwest, suggesting that storage is still a critical component in determining hydrologic regime. Both removal and addition of watershed storage elements is affected by human activities. Thus, it is possible that not all combinations of land-use and watershed storage classes will occur within a region, as was found for the Western Allegheny Plateau ecoregion. If examination of all potential interactions between hydrologic and land-use classes is necessary or desired, we will need to include data from adjacent states with a wider gradient of watershed conditions. The use of nonlinear regression analysis followed by visual analysis of bi-plots to identify variables of interest for hydrologic thresholds is less objective and will be of limited use in cases of collinearity of independent variables, or where classes have alternative definitions, e.g., Condition B < X OR Condition C < Y, as was observed for the Western Allegheny Plateau ecoregion. In these 50 ------- cases, CART analysis should be a more powerful approach to define thresholds. It is possible, however, that graphical analysis will be more useful when interactive terms (Condition B < X and Condition C < Y) are important in denning thresholds. In the current analysis, visual graphical analysis was used because of the form of the equation relating peak flows and watershed variables; no inflection point can be defined in an exponential curve because the second derivative is constant. Graphical analysis can be followed by CART analysis to confirm the placement of the threshold in an unbiased fashion. ACKNOWLEDGMENTS We gratefully acknowledge contributions by Federal and contract staff of the USGS EROS Data Center (Sue Greenlee, Kris Verdin, Jay Kost, Dean Tyler) under interagency agreement DW- 14938973 in developing Stage II and Stage III tools used for catchment aggregation and creation of hydrologicaUy-corrected digital elevation models, without which this project would not have been possible. We also gratefully acknowledge the contributions of Sharon Batterman (US EPA MED- Duluth, IAG Project Officer) and GIS staff of US EPA Region III, in particular, Don Evans and Carmen Constantine, in helping to create Stage II and Stage III 10- and 12-digit HUC boundaries for the state of West Virginia; and of Wendy Bkke-Coleman, US EPA Office of Environmental Information, for facilitating this arrangement. We thank Hope Childers, contract staff at US EPA Region III, Wheeling, WV office for providing the mining-modified NLCD land-use coverage for West Virginia. We thank US EPA staff at Region III- Philadelphia, PA and Fort Meade, MD (Tom DeMoss, Tom Pheiffer) and Wheeling, WV offices (Maggie Passmore, Frank Borsuk); state of WV DNR and DEP staff (Dan Cincotta, WV DNR; Pat Campbell, Dave Montali, John Wilts, WV DEP); and staff of the Canaan Valley Institute (Tom DeMoss, Ron Preston) for consultation on all 51 ------- aspects of this project We thank Debra Taylor at the US EPA MED-Duluth for help in data entry and USGS gaging station watershed delineations. We gratefully acknowledge the GIS support provided by OAO Corporation under US EPA FAIR Contract 68-W5-0065, Delivery Order #24; and Computer Sciences Corporation under US EPA FAIR II Contract 68-W01/W02-032, Task Order #024 in watershed characterization and map production. We also thank Peg Pelletier and Jonathan Smith for providing helpful review comments on an earlier draft of this report. ------- REFERENCES Boward, D.M., P.P. Kazyak, S.A. Stranko, M.K. Kurd, and T.P. Prochaska. 1999. From the Mountains to the Sea: The State of Maryland's Freshwater Streams. EPA/903/R-99/023. Maryland Department of Natural Resources, Monitoring and Nontidal Assessment Division, Annapolis, MD. Brazner, J.C., D.K. Tanner, N.E. Detenbeck, S.L. Batterman, S.K. Stark, L.A. Jagger, and V.M. Snarski. 2003. Landscape influences on fish assembkge structure and function in western Lake Superior tributaries. Submitted to Environmental Management. Cincotta, D. 2000. A Small Watershed Characterization, Classification, and Assessment for West Virginia Utilizing EMAP Design and Tools. R-EMAP proposal submitted to US EPA Mid-Continent Ecology Division, Duluth, MN and US EPA Region III, Philadelphia, PA. Daly, C., R.P. Neilson, and D.L. Phillips. 1994. A statistical-topographic model for mapping climatological precipitation over mountainous terrain. Journal of Applied Meteorology, 33:140-158. Detenbeck N.E., D. Cincotta, J.M. Denver, S.K. Greenlee, and A.R. Olsen. 2003a. Watershed- based survey designs. Submitted to Environmental Monitoring and Assessment (accepted with minor modifications). Detenbeck N.E., C.M. Elonen, D.L. Taylor, L.E. Anderson, T.M. Jicha, and S.L. Batterman. 2003b. Effects of hydrogeomorphic region, watershed storage and mature forest on baseflow and snowmelt stream water quality in second-order Lake Superior Basin tributaries. Freshwater Biology 48(5):911-27. Detenbeck, N.E., S.L. Batterman, V.J. Brady, J.C. Brazner, V.M. Snarski, D.L. Taylor, J.A. Thompson, and J.W. Arthur. 2000. A test of watershed classification systems for ecological risk assessment. Environmental Toxicology and Chemistry 19(4):1174-1181. FGDC. 2002. Federal Standards for Delineation of Hydrologic Unit Boundaries. Federal Geographic Data Committee. http://www.ftw.nrcs.usda.gov/huc_data.html Franken, S.K., D.J.Tyler, and K.L. Verdin. 2001. Development of a National Seamless Database of Topography and Hydrologic Derivatives, Paper 730 in Proceedings, 2001 ESRI International User's Conference, San Diego, CA. Frye, P.M. and G.S. Runner. 1970. A Proposed Streamflow Data Program for West Virginia, U.S. Dept. of the Interior, Geological Survey Water Resources Division, Open-file report, Charleston, WV. Heiskary, S.A. and C.B. Wilson 1990. Minnesota lake water quality assessment report. Minnesota Pollution Control Agency, St. Paul, MN. 53 ------- Herlihy A.T., D.P. Larsen, S.G. Paulsen, N.S. Urquhart, and BJ. Rosenbaum. 2000. Designing a spatially balanced, randomized site selection process for regional stream surveys: The EMAP Mid-Atlantic pilot study. "Environmental Monitoring and Assessment, 63:95-113. Jennings, M.E., W.O. Thomas, Jr., and H.C. Riggs. 1993. Nationwide Summary of U.S. Geological Survey Regional Regression Equations for Estimating Magnitude and Frequency of Floods for Ungaged Sites, 1993. U.S Geological Survey, Reston,VA. Water Resources Investigations Report WRI 94-4002. Johnston, C.A., N.E. Detenbeck, and GJ. NiemL 1990. Cumulative impacts of wetland loss on stream water quality and quantity. Biogeochemistry 10:105-141. Legleiter, KJ. 2001. Interagency Development of National Watershed and Subwatershed Hydrologic Units, Paper 492 in Proceedings, 2001 ESRI International User's Conference, San Diego, CA. NVPDC (Northern Virginia Planning District Commission). 1980. In Schueler, R.T. 1984. The importance of imperviousness. Watershed Protection Techniques 1(3):100-111. Omemik, J.M. and A.L. Gallant. 1988. Ecoregions of the upper Midwest states. EPA/600/3- 88/037, U. S. Environmental Protection Agency, Environmental Research Laboratory, Corvallis, OR. Pan, Y., R.J. Stevenson, B.H. Hill, and A.T. Herlihy. 2000. Ecoregions and benthic diatom assembkges in Mid-Atlantic Highlands streams, USA. Journal of the North American Benthological Society, 19:518-540. Richards, R.P. 1990. Measures of flow variability and a new flow-based classification of great lakes tributaries. J. Great Lakes Res. 16:53-70. Rosgen, D. 1996. Applied river morphology. Wildland Hydrology, Pagosa Springs, CO. Runner, G.S. 1980. Hydrologic Data for Runoff Studies on Small Drainage Areas, West Virginia Department of Highways, Research Project 16, Open-File Report 80-560, Charleston, WV. SAS. 1990. SAS/STAT Users Guide. Version 6, Fourth edition. SAS Institute, Gary, NC. Schueler, T.R. 1994. The importance of imperviousness. Watershed Protection Techniques 1 (3):100-111. Scott, A.G. 1984. Analysis of characteristics of simulated flows from small surface-mined and undisturbed Appalachian watersheds in the Tug Fork basin of Kentucky, Virginia, and West Virginia. U.S. Geological Survey, Reston, VA. Water Resources Investigations Report WRI 84-4151. 54 ------- Tasker, G.D., and J.R. Stedinger. 1989. An operational GLS model for hydrologic regression. Journal of Hydrology 111:361-375. US EPA. 2001. The National Costs of the Total Maximum Daily Load Program (Draft Report). U.S. Environmental Protection Agency, Office of Water, Washington, DC EPA/841/D- 01/003. US EPA. 2002. Consolidated Assessment and Listing Methodology: Toward a Compendium of Best Practices. First Edition. U.S. Environmental Protection Agency, Office of Wetlands, Oceans, and Watersheds. Washington, DC (http://www.epa.gov/owow/monitoring/calm.html). US EPA. 2003. National 303(d) List Fact Sheet. (http: / /oaspub.epa.gov/waters/national_rept. control). Verdin, K.L., and J.P. Verdin. 1999. A topological system for delineation and codification of the earth's river basins. /. Hydrol. 218:1-12. Verry, E.S. 1986. Forest harvesting and water: The Lake States experience. Water Resources Bulletin 22:1039-1047. Waite, I.R., A.T. Herlihy, D.P. Larsen, and D.J. Klemm. 2000. Comparing strengths of geographic and nongeographic classifications of stream benthic macroinvertebrates in the Mid-Atlantic Ffcghlands, USA. Journal of the North American Onthological Society, 19:429^-41. Wilkinson, L. 1999. Systat 9.0 Statistics I. SPSS Inc. Chicago, IL. 55 ------- SEPA United States Environmental Protection Agency Please make all necessary changes on the below label, detach or copy, and return to the address in the upper left-hand comer. If you do not wish to receive these reports CHECK HEREfJ ; detach, or copy this cover, and return to the address in the upper left-hand corner. PRESORTED STANDARD POSTAGE & FEES PAID EPA PERMIT No. G-35 ------- |