Potential contaminant sources and other landscape variables summarized for NHDPlus Version 2.1 catchments within the Chesapeake Bay Watershed

Unless otherwise noted, all analyses were performed in ArcMap version 10.3/4, Excel 2010, or R version 3.3.2.

For variables EPA_septic_total, EPA_septic_active, EPA_septic_abd and EPA_septic_pl_abd, septic system locations were provided by the US Environmental Protection Agency. Counts of septic systems were summarized by NHDPlus Version 2.1 catchment by: 1) associating septic facility locations to NHDPlus catchments, and 2) summing the number of each EPA reported septic system type and total septic systems per catchment.

For the CSO_count variable, Combined Sewer Overflow (CSO) locations for Maryland were downloaded from the Maryland.gov Reported Sewer Overflow Database, for New York from Data.gov, and were provided for EPA Region 3 by the National Pollution Discharge Elimination System (NPDES) Permits Branch Water Protection Division. Data from Maryland were edited to exclude records for Sanitary Sewer Overflow and Bypass categories thereby leaving only CSO reported overflows, which were furthered edited to remove duplicate location records and any sites which fell outside the Maryland state boundary. Data points from New York were edited to remove any locations that fell outside the state boundary or which had no recorded latitude and longitude. CSO locations provided by EPA were edited to remove any VOID or INACTIVE status CSOs. A count of CSOs per NHDPlus Version 2.1 catchment were created by: 1) combining all points into one dataset, 2) associating the CSO locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 3) summing the number of total CSOs per NHDPlus Version 2.1 catchment.

For the FRS_Count variable, Facility Registry Services (FRS) records were downloaded in December, 2016 from the EPA Envirofacts site by state for Virginia, Washington DC, Delaware, New York, Pennsylvania, West Virginia, and Maryland. The data reflects up-to-date facility information as of November 2016. A total count of FRS facilities per NHDPlus Version 2.1 catchment was created by: 1) merging all facility locations into one dataset, 2) removing facilities outside of a 50 meter buffer around the Chesapeake Bay Watershed boundary, 3) removing duplicate records using the “RegistryID”, “lat”, and “long” fields, 5) associating the facility locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 6) summing the number of FRS facilities per NHDPlus Version 2.1 catchment.

For the TRI_Count variable, Toxic Release Inventory (TRI) records were downloaded in December, 2016 from the EPA TRI Factsheet site by state for Virginia, Washington DC, Delaware, New York, Pennsylvania, West Virginia, and Maryland. The data represents 2015 facility information from the 2015 National Analysis Dataset (released October 2016). A total count of TRI facilities per NHDPlus Version 2.1 catchment was created by: 1) merging all facility locations into one dataset, 2) removing facilities outside of a 50 meter buffer around the Chesapeake Bay Watershed boundary, 3) removing duplicate records using the “FRS ID”, “lat”, and “long” fields, 5) associating the facility locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 6) summing the number of total TRI facilities per NHDPlus Version 2.1 catchment.

For the DMR_All, DMR_fed, DMR_POTW, and DMR_nonPOTW variables, EPA Discharge Monitoring Report (DMR) data was downloaded per HUC 4 catchment as 4-level comma separated files (csvs) from the EPA DMR Pollutant Loading Tool. Counts of DMR outfalls were summarized by NHDPlus Version 2.1 catchment by: 1) joining all HUC 4 csvs into one dataset, 2) adjusting positive longitude values to negative, 3) removing duplicate records using the "NPDES number", "Facility Name", "Facility Type", "OutfallNumber", "Monitoring lat", and "Monitoring Long" fields, 4) removing outfalls outside of a 50 meter buffer around the Chesapeake Bay Watershed boundary, 5) manually adjusting the latitude and longitude of 19 outfall locations; the latitude and longitude of these outfalls were located in the Chesapeake Bay waters just outside of the Chesapeake Bay NHDPlus version 2.1 catchments in Maryland and Virginia, but when compared to aerial imagery (ESRI Basemap), appeared to be a facility which extended beyond the edge of the NHDPlus Version 2.1 catchment. The location of the outfall was manually adjusted so that the outfall location was properly recorded in the appropriate NHDPlus Version 2.1 catchment, 6) associating the outfall locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 7) summing the number of each DMR outfall type per NHDPlus Version 2.1 catchment.

For the Septic_num variable, the number of septic tanks per NHDPlus Version 2.1 catchment were generated from two datasets provided by Peter Claggett, USGS Eastern Geographic Science Center: 1) a 30x30m raster of Total Housing Units (THU) 2010, and 2) a 30x30m raster of sewer treatment areas for the Chesapeake Bay Watershed. American Community Survey (ACS) Fact Finder was used to download the 2014 five-year all housing data, Table B25024 for census block groups in Virginia, Washington DC, Delaware, New York, Pennsylvania, West Virginia, and Maryland. Tigerline census blocks for each state were also downloaded from the 2014 Tigerline FTP site. Counts of septic tanks were summarized by NHDPlus Version 2.1 catchment by: 1) removing the sewered areas from the total housing units raster, 2) joining the ACS census data to the Tigerline census block shapefiles via the common “GEOID” column, 3) merging all of the states’ census blocks, 4) clipping the census blocks to a 50 meter buffer around the Chesapeake Bay Watershed, 5) calculating the ratio of detached to total houses by dividing the “HD01_VD02” field by the “HD01_VD01” field, and 6) multiplying the total housing units raster by the calculated ratio. Resulting values were visually compared against actual house footprints on the ESRI World Imagery base map in ArcGIS. From visual inspection it became clear that apartment complexes and townhouse rows on septic systems had no more than 2-3 houses max (aka septic systems) per 30m cell. Therefore all cells in the resulting septic raster with a value greater than 2 were removed. The resulting raster was used to sum the total number of septic tanks per NHDPlus Version 2.1 catchment. Final septic counts per catchment were then rounded to the nearest whole number to eliminate fractional counts.

For variables CAFO_COW_TOT, CAFO_SWINE_TOT, CAFO_POULT_TOT, CAFO_OTH_TOT, and CAFO_ALL_TOT, Confined Animal Feeding Operation (CAFO) animal headcounts were obtained from 1) Mark Zolandz, U.S. Environmental Protection Agency Region 3 National Pollution Discharge Elimination System (NPDES) Enforcement Branch, Water Protection Division, 2) Betsy Bowles, Animal Feeding Operations Program Coordinator, Virginia Department of Environmental Quality, 3) Douglas Ashline, Division of Water, New York State Department of Environmental Conservation, and 4) Kristina Heinemann, Agricultural and Decentralized Wastewater Management Coordinator, U.S. Environmental Protection Agency Region 2 Watershed Management Branch. CAFO animal headcounts per animal type and across all animals were summarized by NHDPlus Version 2.1 catchment by: 1) sorting the CAFO locations by operation type, then by animal type with the following final categories: Cow_total (includes beef, dairy, steer, cattle, heifer, and veal), Swine_total (includes swine, pig, hog, finishing swine, gilt, boar, sow, nursing pig, and farrow/feeder), Poultry_total (includes chicken, layer, broiler, hen, turkey, duck, and pullets), and Other (includes horses, sheep, goats, bison, ponies, and mini horse), 2) associating the sorted CAFO locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 3) summing the head counts of total CAFO animals and CAFO animals by type for each NHDPlus Version 2.1 catchment. For catchments with no reported CAFO animal headcounts, a value of 'ND' (no data) was assigned. ND values may be converted to 0's for statistical analyses at the data user's discretion.

For the variables bioskg_2003n, bioskg_2003p, bioskg_2004n, bioskg_2004p, bioskg_2005n, bioskg_2005p, bioskg_2006n, bioskg_2006p, bioskg_2007n, bioskg_2007p, bioskg_2008n, bioskg_2008p, bioskg_2009n, bioskg_2009p, bioskg_2010n, bioskg_2010p, bioskg_2011n, bioskg_2011p, bioskg_2012n, bioskg_2012p, bioskg_2013n, and bioskg_2013p, nitrogen and phosphorus from biosolids applications per county were downloaded for 2003-2013 from the Chesapeake Bay Midpoint Assessment. Land cover from the National Land Cover Dataset (NLCD) was downloaded for 2001, 2006, and 2011 from the Multi-Resolution Land Characteristics Consortium. Biosolids application data was applied to each NLCD year as follows: NLCD 2001 used with 2003 biosolids data; NLCD 2006 used with biosolids data for 2004-2008; NLCD 2011 used with biosolids data from 2009-2013. Nitrogen and phosphorus applications from biosolids in kilograms were summed by NHDPlus Version 2.1 catchment by 1) identifying al the counties that intersect the Chesapeake Bay Watershed boundary and clipping each NLCD dataset to the county boundaries, 2) reclassifying the land cover to cells which represent Pasture/Hay (NLCD value 81) and Row Crop (NLCD value 82) land uses, 3) obtaining a count of the number of cells in the Pasture/Hay and Row Crop categories per county, 4) converting the polygons of county boundaries to rasters, 5) joining the reclassified raster (only Pasture/Hay and Row Crops) raster to the county raster so that the cells identified as Pasture/Hay and Row Crops also contained the “FIPS” code for each county, 6) joining the county biosolids data to Pasture/Hay and Row Crops raster using the “FIPS” code, 7) calculating pounds of nitrogen and phosphorus from biosolids by dividing reported county applications by the count of Pasture/Hay and Row Crop cells per county, 8) creating a new raster for each field (nitrogen or phosphorus from biosolids for the years 2003-2013), 9) summing the pounds of nitrogen/phosphorus applications from the Pasture/Hay and Row Crop cells in each NHDPlus Version 2.1 catchment, and 10) multiplying the final pounds applied per catchment value by 0.453592 to convert pounds to kilograms.

For the variables CDL2002_phyto, CDL2008_phyto, CDL2009_phyto, CDL2010_phyto, CDL2011_phyto, CDL2012_phyto, CDL2013_phyto, CDL2014_phyto, CDL2015_phyto and CDL2016_phyto the U.S. Department of Agriculture Cropland Data Layer (CDL) was downloaded for 2002 and 2008-2015. The percent of land per NHDPlus Version 2.1 catchment covered by phytoestrogenic crops was determined by 1) clipping the CDL layers to the boundary of counties which intersect the Chesapeake Bay Watershed boundary and resampling 2008 and 2009 layers to 30x30 meter cells, 2) reclassifying the values of crop cover types to include only plants identified as having high levels of phytoestrogens (alfalfa, barley, clover/wildflowers, flaxseed, grapes, oats, peanuts, rye, soybeans, hops and wheat); the reclassification values were 2 for crops identified as year-round crop cover and 1 for plants identified as rotating or half-year crop cover, 3) summing the values of the reclassified raster for each NHDPlus Version 2.1 catchment, 4) separately summing the total number of CDL cells from each raster layer per NHDPlus Version 2.1 catchment, and 5) calculating the percent coverage by dividing the sum of the reclassified values (1 or 2 for high phytoestrogen crops) by 2 so that phytoestrogen cover for all year equaled 1, and half year cover equaled 0.5. These reclassified values were then divided by the total number of CDL cells per catchment and multiplied by 100 to achieve percent coverage of phytoestrogenic plants per NHDPlus Version 2.1 catchment.

For the variables hg2003, hg2004, hg2005, hg2006, hg2007, hg2008, hg2009, hg2010, hg2011, hg2012, hg2013, and hg2014, mercury deposition grids were downloaded from the National Atmospheric Deposition Network for 2003-2014. Total micrograms of mercury deposited from the atmosphere per NHDPlus Version 2.1 catchment were calculated by 1) converting each deposition grid from float to integer-type by multiplying by 10,000, 2) converting the rasters to polygons, 3) calculating the percent of each deposition polygon that covers each NHDPlus Version 2.1 catchment, 4) using the percent of grid coverage to calculate the amount of deposition per catchment as deposition (ug/m2/yr) = (deposition/10,000) * percent coverage per NHDPlus Version 2.1 catchment, 5) adding the deposition per coverage area to find the total deposition per NHDPlus Version 2.1 catchment, and 6) multiplying the deposition (ug/m2/yr/catchment) by the NHDPlus Version 2.1 reported catchment area to find total deposition per catchment per year.

For the Landfill_Count variable, New York and Maryland landfill locations were obtained from Peter Claggett, USGS Eastern Geographic Science Center. The EPA Greenhouse Gas Reporting Program (GHGRP) points for 2015 were downloaded from EPA’s FLIGHT Tool using the search parameters: keyword=landfill; year=2015; GHGs=ALL; data type=All Emitters. The data was reported to EPA by facilities as of 08/13/2016. New York locations were downloaded from the New York State open data site, from which a point shapefile of landfill locations was created using the column “East Coordinate” as the X field and “North Coordinate” as the Y field in UTM Zone 18N. Maryland landfill locations were obtained from the 2015 list of Permitted Solid Waste Acceptance Facilities published by the Maryland Department of the Environment. Pennsylvania municipal waste operations were downloaded from the Pennsylvania Spatial Data Access site on 3/6/2017. A count of landfills for each NHDPlus Version 2.1 catchment were generated by 1) cleaning the Pennsylvania municipal waste operations to remove all facility types aside from “Landfill” and “Landfill-Abandoned” facilities, 2) merging the landfill locations from New York, Maryland, EPA, and Pennsylvania into one spatial dataset, 3) associating the facility locations to NHDPlus Version 2.1 catchments in the Chesapeake Bay Watershed, and 6) summing the number of total landfills per NHDPlus Version 2.1 catchment.

For the variables CATN_2002, CATP_2002, CONN_2002, CONP_2002, HOGN_2002, HOGP_2002, OTHN_2002, OTHP_2002, POULTN_2002, POULTP_2002, TOTN_2002, TOTP_2002, UNCONN_2002, UNCONP_2002, CATN_2007, CATP_2007, HOGN_2007, HOGP_2007, OTHN_2007, OTHP_2007, POULTN_2007, POULTP_2007, TOTN_2007, TOTP_2007, CATN_2012, CATP_2012, HOGN_2012, HOGP_2012, OTHN_2012, OTHP_2012, POULTN_2012, POULTP_2012, TOTN_2012, and TOTP_2012, manure data was obtained from Andy Sekellick, USGS. County estimates of nitrogen and phosphorus in manure were made from the USDA Census of Agriculture livestock population category as kilograms of nutrients in manure per year, calculated as the sum of the animal population multiplied by the N or P manure factor multiplied by the life cycle. For data used and more details on the process to determine county-level estimates, see Mueller and Gronberg (2013) and Gronberg and Arnold (2017). Nitrogen and phosphorus from animal manure in kilograms were summed by NHDPlus Version 2.1 catchment by 1) identifying the counties that intersect the Chesapeake Bay Watershed boundary and clipping the NLCD 2001, 2006, and 2011 dataset to the county boundaries; NLCD 2001 was used for 2002 manure data, NLCD 2006 was used for 2007 manure data, and NLCD 2011 was used for 2012 manure data, 2) reclassifying the land cover to cells which represent Grassland/Herbaceous (NLCD value 71) and Pasture/Hay (NLCD value 81), 3) obtaining a count of the number of cells in the Pasture/Hay and Grassland/Herbaceous categories per county, 4) converting the polygons of county boundaries to rasters, 5) joining the reclassified raster (only Pasture/Hay and Grassland/Herbaceous) raster to the county raster so that the cells identified as Pasture/Hay and Grassland/Herbaceous also contained the FIPS code for each county, 6) joining the county manure data to Pasture/Hay and Grassland/Herbaceous raster using the “FIPS” code, 7) calculating kilograms of nitrogen and phosphorus from manure per cell as nitrogen or phosphorus per county divided by the count of Pasture/Hay and Grassland/Herbaceous cells per county, 8) creating a new raster for each field (nitrogen or phosphorus from manure by animal type for the years 2002, 2007, and 2012), and 9) summing the kilograms of nitrogen/phosphorus in manure from the Pasture/Hay and Grassland/Herbaceous cells in each NHDPlus Version 2.1 catchment.

For the variables PCTbarren01, PCTcult01, PCTdev01, PCTfor01, PCTherb01, PCTshrub01, PCTwater01, PCTwet01, PCTbarren06, PCTcult06, PCTdev06, PCTfor06, PCTherb06, PCTshrub06, PCTwater06, PCTwet06, PCTbarren11, PCTcult11, PCTdev11, PCTfor11, PCTherb11, PCTshrub11, PCTwater11, and PCTwet11, the National Land Cover Database (NLCD) for 2001, 2006, and 2011 were downloaded from the Multi-Resolution Land Characteristics Consortium website. Percent of land cover types by NHDPlus Version 2.1 catchment were calculated by 1) clipping the NLCD 2001, 2006, and 2011 datasets to the boundary of counties which intersect the Chesapeake Bay Watershed boundary, 2) reclassifying each raster into separate masks of the following categories: Water (11), Developed (21,22,23,24), Barren (31), Forest (41,42,43), Shrubland (52), Herbaceous (71), Cultivated (81,82), and Wetland (90,95), 3) summing the number of cells in each category, and 4) dividing the number of cells of each land cover category per catchment by the total number of cells per catchment and multiplying by 100 to get percent land cover type per NHDPlus Version 2.1 catchment.

For the variables Barren01_11, Cult01_11, Dev01_11, For01_11, Herb01_11, Shrub01_11, Water01_11, and Wet01_11, percent land cover change for each category was calculated by subtracting the 2001 NLCD percent coverage classes (PCTbarren01, PCTcult01, PCTdev01, PCTfor01, PCTherb01, PCTshrub01, PCTwater01, PCTwet01) from the 2011 NLCD percent coverage classes (PCTbarren11, PCTcult11, PCTdev11,PCTfor11, PCTherb11, PCTshrub11, PCTwater11, PCTwet11) to get a positive or negative percent change from 2001 to 2011 per each NHDPlus Version 2.1 catchment.

For the Tchg01_11 variable, the total percent of land cover change per NHDPlus Version 2.1 catchment from 2001 to 2011 was calculated by 1) subtracting the value of the NLCD 01 cells from the value of the NLCD 2011 cells, 2) reclassifying all cells with a positive or negative value (indicating change) to a value of 1 and 3) summing the total number of changed cells per NHDPlus Version 2.1 catchment, and 4) dividing the number of changed cells by the total number of cells per NHDPlus Version 2.1 catchment to generate a total percent of land cover change per NHDPlus Version 2.1 catchment.

For each year of the variables starting with 24d, atrazine, bif, fipronil, glyphosate, imidacloprid, metalaxyl, Metolachlor, Simazine, and Thiamethoxam, county-level pesticide application estimates were downloaded as tabular data from the National Water-Quality Assessment Program (NAWQA) Pesticide National Synthesis Project. Following methods detailed by Nancy Baker (NAWQA National Pesticide Synthesis Program, personal communication), and the advice of Dean Hively (USGS Eastern Geographic Science Center), the pesticides of interest were split into two categories: crop-specific use pesticides (Atrazine, Bifenthrin, Clothianidin, Fipronil, Imidacloprid, Metalaxyl, Metolachlor, Simazine, Thiamethoxam), and crop and pasture use pesticides (2,4-D and Glyphosate). The high values for estimated pesticide use were then used for the summaries, and pesticide applications (in kilograms) for each NHDPlus Version 2.1 catchment were generated by: 1) identifying all the counties that intersect the Chesapeake Bay Watershed boundary and clipping the NLCD 2001, 2006, and 2011 dataset to the county boundaries, 2) reclassifying the NLCD 2001, 2006, and 2011 raster to Pasture/Hay (NLCD value 81) and Row Crop (NLCD value 82) cells for use with the crop and pasture use pesticides, 3) reclassifying the NLCD 2001, 2006, and 2011 raster to only Row Crop (NLCD value 82) cells for use with the crop-specific use pesticides, 4) summing the number of Pasture/Hay and Row Crop cells per year per NHDPlus Version 2.1 catchment, 5) summing the number of just Pasture/Hay cells per year per NHDPlus Version 2.1 catchment, 6) combining the Pasture/Hay and Row Crop cells to the raster of counties so that the Pasture/Hay and Row Crop have the “FIPS” code (same was done for just the Pasture/Hay cells), 7) joining sums of the number of cells of Pasture/Hay and Row Crop, and just Pasture/Hay to county polygons along with pesticide application data for the application years associated with the NLCD data year (NLCD 2011 used with pesticides for 2009-2015, NLCD 2006 used for pesticides for 2004-2008, and NLCD 2001 used for pesticides for 2001 and 2003), 8) dividing the county pesticide application by the number of cells of Pasture/Hay and Row Crops or just Pasture/Hay per county, 9) creating a separate raster for each pesticide for each year of reported applications, and 10) summing pesticide application by year for each NHDPlus Version 2.1 catchment to find total kilograms applied per catchment.