U.S. Department of the Interior

U.S. Geological Survey

Scientific Investigations Report 2013–5157

Prepared in cooperation with the National Park Service

Synthesis and Interpretation of Surface-Water Quality and 

Aquatic Biota Data Collected in Shenandoah National Park, 

Virginia, 1979–2009
Cover. Falls on Hogcamp Branch, Shenandoah National Park.
Synthesis and Interpretation of Surface-Water 

Quality and Aquatic Biota Data Collected in 

Shenandoah National Park, Virginia, 1979–2009

By John D. Jastram, Craig D. Snyder, Nathaniel P. Hitt, and Karen C. Rice 

Prepared in cooperation with the National Park Service

Scientific Investigations Report 2013–5157

U.S. Department of the Interior

U.S. Geological Survey
U.S. Department of the Interior

SALLY JEWELL, Secretary

U.S. Geological Survey

Suzette M. Kimball, Acting Director

U.S. Geological Survey, Reston, Virginia: 2013

For more information on the USGS—the Federal source for science about the Earth, its natural and living 

resources, natural hazards, and the environment, visit http://www.usgs.gov or call 1–888–ASK–USGS.

For an overview of USGS information products, including maps, imagery, and publications, 

visit http://www.usgs.gov/pubprod

To order this and other USGS information products, visit http://store.usgs.gov

Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the 

U.S. Government.

Although this report is in the public domain, permission must be secured from the individual copyright owners to 

reproduce any copyrighted materials contained within this report.

Suggested citation:

Jastram, J.D., Snyder, C.D., Hitt, N.P., and Rice, K.C., 2013, Synthesis and interpretation of surface-water quality and 

aquatic biota data collected in Shenandoah National Park, Virginia, 1979–2009: U.S. Geological Survey Scientific 

Investigations Report 2013–5157, 77 p., http://pubs.usgs.gov/sir/2013/5157/.
iii

Acknowledgments

Jeb Wofford, National Park Service, provided detailed information on history, methods, and 

objectives of the aquatic macroinvertebrate and fish monitoring efforts in Shenandoah National 

Park, and his extensive feedback on the analyses presented resulted in numerous improvements 

to the report. Alan Williams, National Park Service, provided monitoring data and spatial data 

collected by the park, and was extremely helpful in providing background information on how 

data were historically collected and stored by park staff. 

The compilation of the comprehensive NPSTORET database and the spatial data was performed 

by USGS employee Tom Burley, whose efforts to restructure the datasets in order to incorporate 

unprecedented data volumes into NPSTORET were vital to the success of the data compilation 

component of this study. 

Four thorough reviews were provided by James Webber, Melissa Riskin, Ellen Porter, and 

Tamara Blett, whose comments greatly improved the clarity and continuity of the report.

v

Contents

Acknowledgments........................................................................................................................................iii

Abstract...........................................................................................................................................................1

Introduction.....................................................................................................................................................2

Purpose and Scope ..............................................................................................................................2

Description of Study Area ...................................................................................................................2

Aquatic-Resource Concerns ..............................................................................................................5

Aquatic Monitoring Datasets..............................................................................................................5

Compilation of the Comprehensive Database...........................................................................................6

Compilation and Analysis of Spatial Data..................................................................................................7

General Approach of Analysis of Spatial Patterns and Temporal Trends in 

Water Quality, Aquatic Macroinvertebrates, and Fish ..............................................................8

Analysis of Water-Quality Data ...................................................................................................................9

Water-Quality Data Structure and Analysis.....................................................................................9

Accounting for Streamflow in Water-Quality Trend Assessments ............................................11

Temporal Trends in Water Quality....................................................................................................12

Determination of Appropriate Trend-Model Formulations..................................................14

Temporal Trends in Water-Quality Constituents...................................................................15

Practical Significance of Trend Analyses.....................................................................18

Temporal Trends in Water Temperature ................................................................................19

Comparison of Water-Temperature Trends and Air-Temperature Trends........................23

Spatial Patterns in Water Quality.....................................................................................................24

Analysis of Aquatic Macroinvertebrate Data .........................................................................................27

Aquatic Macroinvertebrate Data Structure and Analysis...........................................................27

Aquatic Macroinvertebrate Community Characteristics.............................................................29

Spatial Patterns in Aquatic Macroinvertebrates ..........................................................................35

Temporal Trends in Aquatic Macroinvertebrates .........................................................................39

Analysis of Fish Data...................................................................................................................................46

Fish Data Structure and Analysis.....................................................................................................46

Fish Community Characteristics.......................................................................................................48

Temporal Trends and Spatial Patterns in Fish Communities .......................................................49

Temporal Trends and Spatial Patterns in Brook Trout Population..............................................54

Relations Among Environmental Factors, Water Quality, Aquatic 

Macroinvertebrates, and Fish .....................................................................................................60

Stream Acidification...........................................................................................................................60

Water Temperature ............................................................................................................................61

Forest Defoliation................................................................................................................................62

Considerations for Future Monitoring ......................................................................................................63

Water-Quality Monitoring Program .................................................................................................64

Aquatic Macroinvertebrate Monitoring Program .........................................................................64

Fish Monitoring Program...................................................................................................................65

Summary........................................................................................................................................................66

References Cited..........................................................................................................................................67

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring 

Program sites between 1990 and 2009, and 54 sites in 1995...................................................73
vi

Figures

1–3. Maps showing—

 1. Location of Shenandoah National Park in Virginia..........................................................3

 2. Water-quality monitoring stations used for Shenandoah Watershed Study, 

Shenandoah National Park Aquatic Monitoring Program, and the 

Springs and Headwater Streams Study, with underlying geology................................4

 3. Monitoring stations used for temporal trend analysis of 

water-quality constituents .................................................................................................10

4–11. Graphs showing—

 4. Cumulative distribution frequency of daily mean streamflow from five 

streamgages used in Shenandoah Watershed Study in 

Shenandoah National Park................................................................................................12

 5. Relations between concentration trend slope and flow-adjusted trend 

slope for acid neutralizing capacity, pH, and sulfate ....................................................15

 6. Trend-slope p-values for 4-season sampling definition against p-values 

for 12-season sampling definition.....................................................................................16

 7. Acid neutralizing capacity trend slopes and significance by geologic 

category and analysis period ............................................................................................17

 8. pH trend slopes and significance by geologic category and analysis period ..........17

 9. Sulfate trend slopes and significance by geologic category and 

analysis period.....................................................................................................................18

 10. Trend slopes and significance of annual water-temperature 

characteristics by analysis period ...................................................................................19

 11. Time series of mean, minimum, and maximum annual water 

temperature for 14 sites, with period of defoliation 

indicated by vertical dotted lines......................................................................................20

12. Map showing U.S. Geological Survey streamgages used in 

water-temperature trend analysis...........................................................................................21

13. Graphs showing annual water-temperature trend slopes and significance for 

nine U.S. Geological Survey streamgages .............................................................................22

14. Graphs showing annual air-temperature characteristics at the Big Meadows 

and Luray meteorological stations, and annual mean water temperature at the 

Shenandoah River at Lynnwood, Virginia, streamgage, with trend lines .........................23

15. Maps showing spatial patterns in pH, acid neutralizing capacity, sulfate, and 

specific conductance, and underlying geology, in Shenandoah National Park..............25

16. Scatterplot of principal component axes against water-quality variables, 

with significant correlations highlighted ................................................................................26

17. Map showing macroinvertebrate sampling sites in Shenandoah National Park............28

18–27. Graphs showing—

 18. Watershed area of macroinvertebrate sampling sites used for 

temporal and spatial analyses, by geologic class .........................................................29

 19. Taxonomic and trophic composition of aquatic macroinvertebrates 

by sample year computed from samples collected from 24 sites ...............................33
vii

 20. Aquatic macroinvertebrate metric values by geologic class for 

12 metrics significantly influenced by geology ..............................................................36

 21. Percent abundance of Leuctra relative to the overall community, the 

Ephemeroptera, Plecoptera, and Trichoptera assemblage, the intolerant 

taxa assemblage, and the shredder assemblage by watershed type........................37

 22. Relations between aquatic macroinvertebrate metrics and watershed 

area for the six metrics significantly influenced by watershed area .........................39

 23. Shredder% computed with Leuctra, and without Leuctra, 

against watershed area, by geologic class ....................................................................40

 24. Trend slopes for the seven aquatic macroinvertebrate metrics with trends 

independent of geology and watershed area (park-wide) by time period ................42

 25. Trend slopes by geologic class for the seven metrics with long-term 

trend slopes indicating significant geologic effects .....................................................43

 26. Trend slope as a function watershed area for selected aquatic 

macroinvertebrate metrics ................................................................................................44

 27. The number of years required to detect a statistically significant change 

in metric means for metrics that showed significant trends park-wide, 

and metrics with trends dependent on geology.............................................................45

28. Map showing all fish sampling site locations, and those selected for brook 

trout analyses and community analyses, within Shenandoah National Park ..................47

29–38. Graphs showing—

 29. The temporal variation in average abundance of species within 

Shenandoah National Park community focal sites........................................................50

 30. The temporal variation in fish species richness for community focal 

sites by National Park Service site code.........................................................................51

 31. Distributions of mean fish richness by geologic class and against 

watershed area for three time periods analyzed...........................................................52

 32. Fish richness trend slope as functions of geologic class and watershed 

area for the three time periods analyzed ........................................................................53

 33. Distributions of brook trout young-of-year mean abundance by geologic 

class and against watershed area for three time periods analyzed ..........................55

 34. Brook trout young-of-year abundance trend slopes by geologic class 

and against watershed area for three time periods analyzed.....................................56

 35. Distributions of brook trout adult mean abundance by geologic class 

and against watershed area for three time periods analyzed.....................................57

 36. Brook trout adult abundance trend slopes by geologic class and 

against watershed area for three time periods analyzed.............................................58

 37. Per capita rate of increase for adult brook trout at 32 focal sites ..............................59

 38. Z-scores of mean annual temperature and abundance of the genus 

Leuctra for sites with records spanning the entire monitoring period, 

grouped by Shenandoah National Park district.............................................................61
viii

Tables

1. Spatial attributes determined for monitoring stations and associated 

watersheds and the sources of those data..............................................................................7

2. Results of principal component analysis on environmental factors for 

monitored watersheds.................................................................................................................9

3. Shenandoah Watershed Study monitoring stations, periods of record, 

years of record, and number of water-quality results contained in the 

dataset used for trend analysis................................................................................................11

4. Watershed area, geology, and period of record for five streamgages in 

Shenandoah National Park used to determine streamflow statistics ...............................12

5. Correlation coefficients for non-exceedance probabilities from the five 

Shenandoah Watershed Study (SWAS) streamgages and the mean 

non-exceedance probability computed from those five streamgages..............................13

6. Statistics for regression models using various flow statistics to 

estimate acid-neutralizing capacity, by site...........................................................................13

7. Watershed area, mean watershed elevation, and period of record for 

U.S. Geological Survey streamgages used for water-temperature trend analysis.........20

8. Years in which samples were collected for 24 Long Term Ecological 

Monitoring sites used in temporal trend analysis.................................................................30

9. Site information, watershed area, predominant geology, and site elevation 

for sites used in analysis of temporal trends and spatial patterns in 

aquatic macroinvertebrates .....................................................................................................31

10. Description of the 14 macroinvertebrate community metrics evaluated..........................32

11. Factor loadings from principal components analysis of aquatic macroinvertebrate 

metrics for data used in the spatial and temporal analyses of aquatic 

macroinvertebrate communities..............................................................................................34

12. Partial coefficients of determination for variables found to be statistically 

significant predictors for each metric, and overall model R2, from 

interactive stepwise general linear modeling used to assess the 

effects of geology and watershed area on spatial variation 

in macroinvertebrate metrics ...................................................................................................35

13. Description of temporal trends of 14 aquatic macroinvertebrate metrics........................40

14. Number of fish community sampling events per month from 1996 to 2009 ......................46

15. Fish species and abundances reported from 1996 to 2009 from 3-pass 

backpack electrofishing surveys at 107 sites ........................................................................48

16. Stepwise general linear model results for effects of geology and watershed 

area on stream fish community richness, brook trout young-of-year abundance, 

and brook trout adult abundance among time periods of 5, 10, and 14 years..................49
ix

Conversion Factors

Multiply By To obtain

Length

centimeter (cm) 0.3937 inch (in.)

micrometer (μm) 0.003937 inch (in.)

meter (m) 3.281 foot (ft) 

kilometer (km) 0.6214 mile (mi)

kilometer (km) 0.5400 mile, nautical (nmi) 

meter (m) 1.094 yard (yd) 

Area

hectare (ha) 2.471 acre

square kilometer (km2) 247.1 acre

hectare (ha) 0.003861 square mile (mi2) 

square kilometer (km2) 0.3861 square mile (mi2)

Flow rate

cubic meter per second (m3/s) 35.31 cubic foot per second (ft3/s)

Temperature in degrees Celsius (°C) may be converted to degrees Fahrenheit (°F) as follows:

°F=(1.8×°C)+32

Horizontal coordinate information is referenced to the North American Datum of 1983 (NAD 83).

Specific conductance is given in microsiemens per centimeter at 25 degrees Celsius (μS/cm at 25 °C).

Concentrations of chemical constituents in water are given in microequivalents per liter (μeq/L).
x

Abbreviations

age 1+ age greater than 1 year

AM aquatic macroinvertebrate

ANC acid neutralizing capacity

CI confidence interval

DOM2 relative abundance (percent) of the two most abundant taxa

DPAC distribute parents among children

E% percent abundance of individuals in order Ephemeroptera

EPT orders Ephemeroptera, Plecoptera, and Trichoptera

EPT% percent abundance of individuals of orders Ephemeroptera, Plecoptera, and 

Trichoptera

EPT-r total number of taxa in orders Ephemeroptera, Plecoptera, and Trichoptera

GLM general linear modeling

Hapto% percent abundance of Haptobenthos (macroinvertebrates requiring clean, firm, 

coarse substrates)

Hydro:T% percent abundance of caddis flies (Trichoptera) in family Hydropsychidae 

(Trichoptera)

IDAS Invertebrate Data Analysis System 

Intol% percent abundance of macroinvertebrates with tolerance values <2 

Leuctra:P% percent abundance of stoneflies (Plecoptera) represented by Leuctra

LOWESS locally weighted scatterplot smoothing

LTEM Long Term Ecological Monitoring (program)

MDC minimum detectable change

NO3

– nitrate

NOx nitrous oxides

NPS National Park Service

PCA principal components analysis

PIBS portable invertebrate box sampler

PTV pollution tolerance value

R2 coefficient of determination

RMSE root mean squared error

SCI stream condition index

Scraper% percent abundance of functional feeding group containing “scrapers” 

Shredders% percent abundance of functional feeding group containing “shredders” 

SHS Springs and Headwater Streams (Study)

SimpsonD Simpson diversity index (1-D)

SO2 sulfur dioxide

SO4

2– sulfate

SWAS Shenandoah Watershed Study

USGS U.S. Geological Survey

YOY young-of-year

YTDC years to detect a significant change
Abstract

Shenandoah National Park in northern and central 

Virginia protects 777 square kilometers of mountain terrain 

in the Blue Ridge physiographic province and more than 

90 streams containing diverse aquatic biota. Park managers 

and visitors are interested in the water quality of park streams 

and its ability to support healthy coldwater communities and 

species, such as the native brook trout (Salvelinus fontinalis), 

that are at risk in the eastern United States. Despite protection 

from local stressors, however, the water quality of streams 

in the park is at risk from many regional stressors, including 

atmospheric pollution, decline in the health of the surrounding forests because of invasive forest pests, and global climate 

change. In 2010, the U.S. Geological Survey, in cooperation 

with the National Park Service, undertook a study to compile, analyze, and synthesize available data on water quality, aquatic macroinvertebrates, and fish within Shenandoah 

National Park. Specifically, the effort focused on creating a 

comprehensive water-resources database for the park that can 

be used to evaluate temporal trends and spatial patterns in the 

available data, and characterizing those data to better understand interrelations among water quality, aquatic macroinvertebrates, fish, and the landscape.

Data from three primary sources, namely the Shenandoah Watershed Study, the Shenandoah National Park Aquatic 

Macroinvertebrate Monitoring Program, and the Springs 

and Headwater Streams Study, were compiled and loaded 

into the National Park Service’s NPSTORET database. This 

effort yielded a comprehensive database containing nearly 

1.3 million measurements of habitat characteristics, approximately 442,000 measurements of water-quality characteristics, 

and over 438,000 measurements of biological taxa (fish and 

aquatic macroinvertebrates), collected across 673 sites over a 

period of more than 30 years.

Temporal trends in water quality indicate conflicting patterns in terms of acidity. Long-term (20- and 30-year) trends in 

acid-neutralizing capacity (ANC) and pH may indicate some 

improvement (decreasing acidity), but short-term (5-year) 

trends suggest increasing acidity. The long-term increase in 

pH occurred park-wide, although the increases were minimal 

in watersheds having siliciclastic bedrock. The long-term 

increases in ANC were mostly limited to watersheds with 

basaltic bedrock. Trends in concentrations of stream-water sulfate, another constituent of atmospheric deposition, indicated 

long-term improvements (declines in concentration) in watersheds having basaltic bedrock, long-term increases in concentration in watersheds having granitic bedrock, but no trend 

in watersheds with siliciclastic bedrock. Park-wide increases 

in mean, median, and maximum water temperatures were 

detected over the last 20 years. The average annual increase 

in mean water temperature park-wide was 0.04 °C, which 

equates to about 1.2 °C over the last 30 years. Short-term 

trends generally coincided with long-term trends but were 

more variable. Water temperatures generally tracked air temperatures, and additional analyses of longer-term (greater than 

80 years) regional air-temperature data showed that the most 

recent increases in air temperature are not unprecedented.

Analysis of spatial patterns in water quality demonstrated 

that watersheds with higher mean elevations, lower land-surface gradients, and greater proportions of basaltic and carbonate geology are least affected by acidification and tend to be 

improving over time. Watersheds having greater proportions of 

siliciclastic and granitic geology, with smaller watershed areas 

and higher minimum watershed elevation tend to be affected 

by acidification and are experiencing continued degradation 

in water quality. There was no apparent spatial pattern in the 

water-temperature trends. 

Benthic macroinvertebrate community metrics were 

found to be highly correlated with geology and, to a lesser 

extent, watershed area. Temporal trends in benthic macroinvertebrates showed evidence of change in community 

structure over time, which in most cases indicated declines 

in stream condition. Although the overall condition of park 

streams would be considered by most measures to be relatively healthy, streams in siliciclastic watersheds, in particular, 

have been and continue to be affected by acidic deposition. In 

addition, park streams have warmed significantly over the last 

20 years and evidence indicates that benthic macroinvertebrate 

communities have responded to the warming trend. 

Analysis of temporal trends in fish-community metrics 

revealed increasing species richness from 1996 to 2010, 

with the greatest increases observed in the larger streams. 

Synthesis and Interpretation of Surface-Water Quality 

and Aquatic Biota Data Collected in Shenandoah 

National Park, Virginia, 1979–2009

By John D. Jastram, Craig D. Snyder, Nathaniel P. Hitt, and Karen C. Rice 
2 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Furthermore, increases in fish richness were detected only 

in watersheds underlain by granitic and basaltic geology; 

no trend in richness was detected in siliciclastic watersheds, 

which also had the lowest richness values. The low richness 

values and lack of improvement over time in siliciclastic 

watersheds further suggests that recovery from acidification 

is not yet occurring in these areas. Increases in richness may 

partially be due to increased water temperatures.

 Analysis of brook-trout population data indicated a 

response to acidification in siliciclastic watersheds, as mean 

abundances of young-of-year (YOY) and adult (age greater 

than 1 year; age 1+) fish were consistently lower in these 

watersheds than in others. Long-term (15-year) increases were 

detected, however, in adult abundances in the siliciclastic 

watersheds, although YOY abundances showed a declining 

trend during this same period. Despite observed demographic 

and environmental variation, brook trout population growth 

rates were stable throughout the study area between 1996 

and 2009.

Although the current monitoring design in the park is 

spatially and temporally robust, specifically in terms of the 

ability to detect changes in water quality and aquatic fauna, 

and to understand processes related to stream acidification, 

multiple potential improvements to the monitoring programs 

were identified. Potential changes that would allow for more 

efficient accomplishment of monitoring objectives or more 

complete representation of water-resources and aquatic fauna 

conditions include reducing the frequency of water-quality 

trend monitoring to quarterly sampling, operating additional 

streamgages within the park, and co-locating monitoring 

stations used for water-quality, aquatic macroinvertebrate, 

and fish monitoring programs. Although these changes would 

potentially reduce monitoring costs or improve the datasets, 

the current dataset was found to be sufficient to satisfy multiple objectives, including objectives or data analysis for which 

the monitoring program was not originally designed, such as 

water-temperature trend analysis. 

Introduction

Shenandoah National Park is located within the Blue 

Ridge physiographic province in northern and central Virginia 

and protects 777 square kilometers (km2) of mountain terrain 

in the Blue Ridge physiographic province. The mountain-top 

location of the park affords protection of more than 90 streams 

containing diverse aquatic biota. Park managers and visitors 

are interested in the water quality of streams in the park, in 

particular, its ability to support healthy populations of native 

brook trout that are at risk in the eastern United States. Despite 

protection from local stressors, however, the water quality of 

streams in the park is at risk from various stressors, including 

atmospheric pollution, declines in the health of the surrounding forests because of invasive forest pests, and global 

climate change.

In 2010, the U.S. Geological Survey (USGS), in cooperation with the National Park Service (NPS), initiated a study to 

examine three separate databases that contain data collected 

within Shenandoah National Park, which contain water-quality 

data, aquatic macroinvertebrate (AM) data, and fish data. 

The purpose of this study was to analyze the three existing 

databases, synthesize the databases, and interpret the data in 

the context of the effect of water quality on the distribution of 

biota in park streams. More specifically, the objectives of this 

study were to:

1. Determine whether relations exist among aquatic faunal distribution and abundance and water quality in park streams 

and assess the relation of those data to ecosystem data;

2. Determine whether the monitoring programs are adequate 

to characterize the status and trends in water quality of 

streams and aquatic biota across the spectrum of acid 

sensitivity within the park; and

3. Analyze and improve the park’s water-resources monitoring database.

Purpose and Scope

The purpose of this report is to describe aquatic resource 

data that have been collected in Shenandoah National Park 

over the last three decades, and to describe a recent effort to 

compile and analyze those data. The data, which include measures of water quality, aquatic macroinvertebrates, and fish, 

were collected by multiple entities to satisfy varied objectives, 

and the analyses reported herein are the result of the first time 

these disparate datasets have been analyzed collectively. This 

report provides an understanding of how water quality, aquatic 

macroinvertebrates, and fish have changed over the period of 

the study and how they vary spatially within the park. Additionally, this report presents analyses of the interrelations 

among the temporal and spatial patterns observed in water 

quality, aquatic macroinvertebrates, and fish.

Description of Study Area

Shenandoah National Park is a narrow strip of NPS land 

169 kilometers (km) in length, located in northern Virginia 

(Va.) (fig. 1). Of the 777 km2 protected within the park boundary, about 40 percent (311 km2

) is designated as wilderness. 

The park is mostly forested, containing a mix of deciduous 

and coniferous trees. The climate is humid temperate, with 

mean annual air temperatures ranging from 9 degrees Celsius 

(°C) at higher elevations to 12 °C at lower elevations, with a 

mean annual rainfall of 100 to 150 centimeters (cm) per year 

(National Park Service, 2012). 

Water resources in the park include over 1,000 km of 

stream channels contained within about 90 perennial streams 

(National Park Service, 1998) that ultimately drain into Chesapeake Bay through three rivers: the Rappahannock River, the 
Introduction 3

James River

River

Rappahannock 

Potomac 

River

Shenandoah 

River

River

S Fork Shenandoah

Chesapeake Bay

N Fork Shenandoah River

WEST VIRGINIA

KENTUCKY

TENNESSEE

OHIO

PENNSYLVANIA

NORTH CAROLINA

VIRGINIA

MARYLAND

DE

NJ

Base from U.S. Geological Survey digital 

line graph, 1:2,000,000, 1987

Figure 1. Map showing location of Shenandoah National Park in Virginia.

84° 83° 82° 81° 80° 79° 78° 77° 76°

41°

40°

39°

38°

37°

EXPLANATION

Shenandoah National Park

Blue Ridge Physiographic 

 Province

Virginia state boundary

0 50 100 KILOMETERS

0 50 100 MILES

VIRGINIA

Location of

study area

Physiography from Fenneman and 

Johnson (1946), 1:7,000,000

Figure 1. Location of Shenandoah National Park in Virginia.

Potomac River, and the James River. These aquatic habitats, 

along with their associated riparian areas, constitute unique 

and important resources in the park, providing irreplaceable 

habitat for many aquatic and terrestrial species. 

Eleven geologic formations have been identified in the 

park (Virginia Division of Mineral Resources, 2003; Southworth and others, 2009). For the purpose of this study, the 

11 formations were grouped into three classes—siliciclastic, 

granitic, and basaltic—based on their ability to neutralize 

acidic inputs. Each of the three classes represent approximately one-third of the area within the park. The siliciclastic 

class is found mainly in the westernmost parts of the park, the 

basaltic class crops out in the central region of the park, and 

the granitic class is present along the eastern margins of the 

park (fig. 2).

The bedrock control on stream-water quality has long 

been recognized in the Blue Ridge physiographic province 

(Bricker and Rice, 1989) and within the Shenandoah National 

park in particular (Cosby and others, 2006; Deviney and 

others, 2006; Rice and others, 2007). In general, watersheds 

developed on siliciclastic bedrock have streams with very 

low ionic-strength waters and little acid-neutralizing capacity (ANC). Watersheds developed on basaltic bedrock have 

streams with the highest ionic-strength and ANC values in 

the park, and watersheds developed on granitic bedrock have 

streams with ionic-strength and ANC values intermediate 

between those developed on siliciclastic and basaltic bedrock. 

Thus, there is a gradient in sensitivity to inputs by acidic deposition across streams in the park that is strongly controlled by 

the type of bedrock that underlies the watershed.
Figure 2. Maps showing water-quality monitoring stations used for A, Shenandoah Watershed Study (SWAS), B, ShenandNational Park Aquatic Monitoring Program (SHEN), and C, the Springs and Headwater Streams Study (SHS), with underlyigeology. Location of study area provided in figure 1.

!

!

!

!

!

!

!

! 

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!

!

!

!

! !

!

!

!

!!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!

!

! !

!!!

! ! !

!!

!!

!

!

!

!

!

!

!

!

!

! !!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!!

!!

! !

!

!

!!

!!

!

!!

! !

!

! !

!

! ! !!

!

!

!

!

!! ! ! !

! !

!

!

!

!

!!

!

!!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!! !

!

!

!

! !

! !

!

!

!

!!!

!

!

!!

!

!

!

!

!

!!

!!

!!!

!!!

!!!

!!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

! !

!

!!!

!!!!!!!!!! !

!

!!!

!

!

!

!

!

!

!

!

!

!

N

!!

!

S SW 

B

Base from U.S. Geological Survey digital line graph

A B C

0 10 20 KILOMETERS0 10 20 MILES

Figure 2. Water-quality monitoring stations used for A, Shenandoah Watershed Study (SWAS), B, Shenandoah National Park Aquatic Monitoring PC, the Springs and Headwater Streams Study (SHS), with underlying geology. Location of study area provided in figure 1.
Introduction 5

Aquatic-Resource Concerns

Pristine water quality is essential to the park’s ecosystems, as well as to the hundreds of thousands of people in the 

surrounding lowlands who rely on water that originates in the 

park for their domestic water supply. The park’s legislative 

history speaks of “beautiful cascading streams” and the critical 

importance of wildlife, recreation, and mountain scenery to 

the establishment of the park (U.S. Department of the Interior, 

1932). The park’s brook trout habitat is widely recognized as 

being among the best in the middle Appalachians, and streams 

in the park collectively provide habitat for at least 38 fish 

species (National Park Service, 2013) and approximately 

240 macroinvertebrate species (National Park Service, 2008). 

As a result, water quality and the effects of degraded water 

quality on the park’s aquatic ecosystems are some of the most 

significant resource-management issues at the park. 

Currently, one of the largest resource-management 

problems that resource managers at the park face is air quality, 

specifically acidic deposition, and its effects on stream-water 

quality in the park. Acidic deposition and its effects on aquatic 

ecosystems have been major scientific and public policy issues 

for nearly 30 years, because of the potential for large-scale and 

long-term changes in water quality and the resulting loss of 

aquatic biota. Previous studies (Cosby and others, 1991, and 

references therein) have noted that streams in large areas of 

the southeastern United States are vulnerable to acidification, 

and the park has received some of the highest levels of acid 

deposition in the country. Acid deposition has declined since 

implementation of the Clean Air Act’s Acid Rain Program in 

1990; however, recovery of many acidified lakes and streams 

is not likely without further reductions in acid deposition 

(Burns and others, 2011). The Shenandoah National Park 

Resources Management Plan (National Park Service, 1998) 

contains 7 project statements on acidic deposition and its 

effects on the environment, and 13 additional project statements concerning other aspects of air-quality management 

(out of 168 total project statements), showing the importance 

of this issue to the park. Characterizing the link between air 

quality and water quality related to acidic deposition, contaminants, soil conditions, and forest health is a high priority for 

research in Shenandoah National Park. 

In addition to acidic deposition, park managers are 

contending with other issues relevant to aquatic resources, 

including managing potential effects of climate change, 

preservation of native brook trout, introduction of exotic fish, 

and forest defoliation by insect pests. Climate change threatens to alter the temperature and precipitation regimes of the 

park, which in turn could disrupt the balance of the cool- and 

cold-water ecosystems in headwater streams of the park. Many 

of the aquatic organisms that inhabit these streams, including the native brook trout, require cool water temperatures to 

complete their life cycles; therefore, warming of park streams 

poses a risk to these unique communities. Warming stream 

temperatures also may increase the abundance of exotic organisms in park streams as water-temperatures increase to levels 

more acceptable to such organisms, creating an additional 

stress for native species. Changes in the precipitation regime 

as a result of climate change could also threaten these ecosystems as disturbances such as floods and droughts may occur 

more frequently, with greater duration and greater intensity.

Forest defoliation by the gypsy moth (Lymantria dispar) 

has been an important resource-management issue in the park 

over the last three decades. The most noticeable effects of 

defoliation occur in the terrestrial ecosystems, but there are 

potential cascading effects on water quality and aquatic biota, 

centered around changes in the aquatic food web. In-stream 

nutrient concentrations may increase, and the type of available 

food sources may change, as uptake of nutrients by trees is 

decreased and inputs to the stream from leaf fall are decreased. 

These changes, potentially exacerbated by increased light 

levels, could possibly shift the trophic composition of the 

AM community from organisms that consume leaf material 

deposited in the stream to those that consume algae growing 

on stream substrates. 

Aquatic Monitoring Datasets

Multiple long-term and (or) spatially extensive aquatic 

data-collection efforts were identified as primary data sources 

for this study. Although this collection of datasets is not inclusive of all aquatic data collected in the park, it is believed to 

represent the most temporally and spatially extensive datasets 

containing data pertinent to the objectives of this study.

Recognizing the potential threat of acidic deposition, the 

NPS initiated the Shenandoah Watershed Study (SWAS) in 

1979 in cooperation with researchers from the University of 

Virginia. The goal of SWAS has been to monitor and assess 

the acid-base chemistry in park streams in order to better 

understand the effects of acidic deposition on stream resources. 

SWAS data include stream discharge, water temperature, and 

several water-quality parameters relevant to acidic deposition (specific conductance, pH, ANC, and the concentrations 

of major cations and anions). Since the inception of SWAS, 

water-quality data have been collected at least quarterly at sites 

in 14 watersheds (fig. 2A), and multiple samples have been collected during individual storm events at sites in 3 watersheds. 

The SWAS data analyzed for this report include records from 

20 of the 29 sites in the complete SWAS dataset, although the 

number of years sampled varies considerably among sites. This 

dataset contains nearly 195,000 individual water-quality results 

collected during nearly 11,000 sampling events, with over 

150,000 of those results collected at 6 of the 20 sites. 

In addition to SWAS data, the NPS has collected waterquality data as part of their biological monitoring efforts. 

These monitoring data were intended to provide context for 

the fisheries and AM data, and include measurements of 

water temperature, dissolved oxygen, specific conductance, 

pH, and total dissolved solids. The dataset provided for this 

study includes over 9,000 water-quality results from 204 sites 

(fig. 2B), dating back to 1986, with the majority of the results 

coming from approximately one-third of the sites.
6 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Incorporation of biological endpoints into the aquatic 

monitoring program began in 1982 with the start of the 

Shenandoah National Park Fisheries Monitoring Program. 

The initial focus of the program was on monitoring the status 

of the brook trout fishery within the park, although the scope 

of the program was expanded in 1995 to include the occurrence and abundance of all fish species and associated fish 

community metrics (J.E.B. Wofford, National Park Service, 

written commun., 2011). In total, the fish data include over 

370,000 records from over 1,600 sampling events across 

211 sites. 

The AM monitoring program for the park was initiated 

in 1984 when the park was designated as a prototype for the 

Long-Term Ecological Monitoring (LTEM) Program of the 

NPS. The objectives of the AM monitoring program are to use 

AM community and habitat metrics to characterize status and 

trends in stream condition. The AM dataset provided for this 

study contains approximately 65,000 results for 346 measures 

(taxa counts) from 737 sampling events across 114 sites.

The Springs and Headwater Streams (SHS) Study was 

conducted by staff from the USGS Leetown Science Center 

from 2007 through 2008. This study was initiated because of 

a lack of water-quality data collected from headwater streams 

and springs in the park. In addition, the Water Resources 

Scoping Report for the park (Vana-Miller and Weeks, 2004) 

suggested that the SWAS program, which was focused on 

sampling at low elevations in the watersheds, needed to consider expansion to sampling water quality higher in the watersheds. Therefore, the goal of this study was to characterize the 

physical, chemical, and biological characteristics in headwater 

streams and springs in the park in order to determine if current monitoring activities adequately characterize the range 

of aquatic habitat in the park. This goal was achieved through 

sampling water quality and AM communities at 23 headwater 

streams, 40 springs, and 9 LTEM stream sites (fig. 2C). Waterquality data collected for this effort include the same list of 

constituents as those collected for the SWAS program. This 

effort resulted in the collection of nearly 2,500 water-quality 

measurements and approximately 5,750 biological measures 

(taxa counts) from 140 sampling efforts across the sites.

Compilation of the Comprehensive 

Database

The comprehensive database, consisting of fish and AM 

monitoring data collected by Shenandoah National Park, 

water-quality data collected by SWAS, and AM and waterquality data collected by the SHS Study, was compiled in the 

NPSTORET framework of the NPS. This framework was 

selected for consistency with NPS protocols for data storage, and to facilitate data upload to the U.S. Environmental 

Protection Agency’s STORET database. Water-quality data 

collected by as part of the SWAS program, as well as water 

data collected in the park for numerous other studies, already 

existed in an NPSTORET database. This existing NPSTORET 

database was used as the foundation to which the other 

datasets were added. 

Water-quality data and site information from the SHS 

study were previously loaded in the NPSTORET database; 

therefore, only the AM data collected as part of that study 

remained to be loaded into the NPSTORET database. Counts 

of AM taxa per sampling event, as provided by the study 

authors, were reformatted into row-major format and loaded 

into NPSTORET.

Water-quality, AM, fish, and habitat data collected by 

the NPS are maintained by the NPS in a Microsoft Access 

database. Prior to loading the data into NPSTORET, these 

data required extensive reformatting to collapse the multiple 

column-major tables, originally linked by one-to-many relationships, into fewer row-major tables, while still maintaining 

the integrity of relationships used to link data to the appropriate sampling event and location. The resultant NPSTORET 

database contained nearly 1.3 million measurements of habitat 

characteristics, approximately 442,000 water-quality measurements, and over 438,000 measurements of biological taxa 

(AM and fish data) collected across 673 sites. This comprehensive database formed the foundation for the analyses presented in this report, and it is expected to become the primary 

database for archival and analysis of all historic and future 

aquatic monitoring data collected in the park.

Concentration units of water-quality constituents were 

standardized for storage in NPSTORET such that concentrations are stored in units of microequivalents per liter (μeq/L). 

Therefore, concentrations of chemical constituents are presented in μeq/L throughout this report.
Compilation and Analysis of Spatial Data 7

Compilation and Analysis of 

Spatial Data

In addition to the water-quality and biological (AM and 

fish) data just described, a suite of other environmental variables were used as covariates in spatial and temporal trend 

analyses in order to determine the extent to which spatial patterns and temporal trends in water quality and biological measures depend on local habitat or landscape setting. The selected 

environmental variables included those known or expected to 

influence water quality and biota, such as geology, elevation, 

and watershed area, as well as variables for which data were 

readily available and that may be of interest in the current study 

or potential future studies (for example, forest vegetation composition, location of roads and trails, and stream order). 

Environmental data were compiled from various sources 

including published datasets specific to the park, published 

datasets with Nation-wide coverage, and datasets maintained 

by park staff (table 1). Local site conditions were characterized 

for all sites, and a subset of the watersheds was also characterized to determine the landscape attributes associated with each 

site. The subset characterized consisted of those watersheds 

that had sites with sufficient monitoring data for trend analysis. 

 The datasets just described were manipulated using 

ArcGIS 9.3.1 to calculate and assign attributes to the sampling points or watershed areas draining sampling points. In 

addition, watershed boundaries were delineated and attributes 

derived for all unmonitored streams leaving the boundaries 

of the park to provide a basis for extrapolating analyses to 

unmonitored watersheds should that objective be part of any 

future work. The resulting collection of spatial data was compiled into an ArcGIS 9.3.1 geodatabase (table 1). 

Table 1. Spatial attributes determined for monitoring stations and associated watersheds and the sources of those data.

Attribute Units Source

Watershed

Area Square kilometers Measured

Elevation

 minimum

 maximum

 mean

Meters National Elevation Dataset

Slope

 mean

Percent National Elevation Dataset

Aspect Percent in direction classification National Elevation Dataset

Land cover Percent per classification 2001 National Land Cover Dataset (Homer and others, 2004)

Vegetation Percent per classification Shenandoah National Park internal data

Gypsy moth defoliation Presence or absence of defoliation by year Shenandoah National Park internal data

Mean precipitation (1980–2009) Millimeters PRISM Climate Group

Roads and trails Percent cover Shenandoah National Park internal data

Bridges and culverts Count within watershed Shenandoah National Park internal data

Road and trail stream crossings Count Shenandoah National Park internal data

Wastewater treatment facilities Presence/absence Shenandoah National Park internal data

Geologic classification Percent Morgan and others (2003)

Geologic units Percent Southworth and others (2009)

Site

Elevation Meters National Elevation Dataset

Aspect Percent in direction classification; 

decimal degrees

National Elevation Dataset

Stream order

 Strahler (1952)

 Shreve (1966)

Order National Hydrography Dataset

Vegetation Percent per classification in 100-meter buffer Shenandoah National Park internal data

Stream slope

 100 meters above site

 NHD reach 

Percent National Hydrography Dataset

Geologic classification Classification Morgan and others (2003)

Geologic units Unit Southworth and others (2009)
8 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

General Approach of Analysis of 

Spatial Patterns and Temporal 

Trends in Water Quality, Aquatic 

Macroinvertebrates, and Fish 

The fish and AM monitoring data collected by park staff 

have not been previously analyzed for long-term trends. In 

addition, there has been no prior synthesis or analysis of these 

data in conjunction with the SWAS water-quality data or the 

SHS Study. The combination of these datasets facilitated an 

integrated assessment of status and trends in aquatic resources 

and provided a basis to generate hypotheses regarding the 

effects of regional anthropogenic stressors on stream condition. 

A consistent analysis approach was used to evaluate 

spatial patterns and temporal trends in water quality, AM, 

and fish. Specific characteristics of the individual datasets 

(water quality, AM, and fish), however, required that different 

statistical techniques be applied to each dataset. The specific 

methods of analysis, therefore, are presented in the respective 

sections of this report to maintain clarity in how the analyses 

were performed.

There were several discrepancies among the three data 

types (water quality, AM, and fish) that prevented a fully integrated data analysis. Specifically, the location of sample sites 

varied among monitoring programs, as did the timing and frequency of sampling, preventing direct assessment of relations 

between different data types. Therefore, the general analytical 

approach taken was to assess spatial patterns and temporal 

trends for each data type individually and subsequently look 

for commonalities in short-term and long-term trends in 

order to derive hypotheses regarding causal explanations for 

observed patterns. In particular, temporal trends were assessed 

in the context of acidic deposition, forest defoliation caused by 

the gypsy moth, and climate change—the three stressors most 

likely to influence stream conditions in the park. 

Numerous response variables were assessed for each 

data type. Water-quality response variables included mean, 

maximum, and minimum values for water temperature, and 

numerous water-quality constituents known to be important 

to drivers of water quality in the park, including pH, ANC, 

sulfate (SO4

2–), and nitrate (NO3–

). Aquatic macroinvertebrate response variables included 15 metrics that describe the 

diversity and structure of AM communities and often are used 

to infer stream condition. Fish response variables included 

3 measures of the abundance and demography of brook trout 

as well as 22 measures of the diversity and composition of the 

overall fish community. 

Spatial patterns and temporal trends were assessed 

through a multiple regression modeling approach that related 

each individual response variable to a suite of potential 

predictor variables. For analysis of spatial patterns, it was not 

feasible to evaluate the relations among the large number of 

site and watershed characteristics and the suite of response 

measures. There were numerous correlations among the lists 

of environmental variables; therefore, principal components 

analysis (PCA; Shaw, 2003) was used to reduce the list of 

predictor variables to a smaller number of uncorrelated multivariate gradients. The PCA results showed that environmental 

factors can be reduced to 3 statistically significant (p<0.10) 

dimensions that collectively describe nearly 62 percent of the 

variability in the environmental characteristics of these watersheds (table 2). The first PCA axis explains over 35 percent 

of the variation and is highly correlated with the proportion 

of siliciclastic geology in the watershed and mean watershed 

slope. The second PCA axis explains an additional 14 percent 

of the variation and is most strongly correlated with watershed 

area and minimum elevation. The third PCA axis explains an 

additional 12 percent of the variation and positively correlates 

with the proportion of granitic geology in the watershed and 

negatively correlates with the proportion of basaltic geology. 

The strong correlations between vegetation characteristics and 

individual principal components is understood to be a result 

of the correlation between those vegetation characteristics and 

other physical attributes (such as geology and elevation), and 

therefore the vegetation characteristics would be expected to 

have associative relations, as opposed to causative relations, 

with water quality, fish, and AM measures. Thus, geology and 

watershed area, and their interaction, were used as predictor 

variables in spatial analyses.

Temporal trend analyses focused on changes in response 

measures through time, using time (in years) as the primary 

predictor variable. Trend analyses were performed for four 

periods; namely, the most recent 20-year period (1990–2009), 

14-year period (1996–2009), 10-year period (2000–2009), 

and 5-year period (2005–2009). These periods were selected 

to represent changes over time that have led to currently 

observed conditions. Selection of different time periods and 

(or) start/end dates could possibly yield different results. The 

analysis of fish data was limited to collections made after 

1995 because of concerns about the comparability of sampling methods over time, hence the 14-year period instead 

of a 15-year period consistent with the 5-year interval of 

the other analyses. Aquatic macroinvertebrate and fish data 

were collected at an annual time step; therefore, simple linear 

regression was used to evaluate changes in response measures 

over time. Monthly water-temperature data were summarized 

as annual statistics, and those annual statistics were analyzed 

similarly to AM and fish data. Water-quality data were collected more frequently and therefore were analyzed using the 

Seasonal Kendall Trend Test (Helsel and Hirsch, 2002). Multiple linear regression was subsequently used to evaluate the 

importance of watershed characteristics (geology, watershed 

area, and the interaction between geology and watershed area) 

on temporal changes in response variables (that is, regression slopes) for each time period. In general, the analytical 

approach just described was used to assess spatial and temporal trends in response to water-quality, AM, and fish variables.
Analysis of Water-Quality Data 9

Analysis of Water-Quality Data

Water-quality data were compiled and analyzed to evaluate 

temporal trends and spatial patterns in water-quality constituents 

related to stream acidification. The temporal trend analysis was 

developed in a manner that provides an understanding of shortterm and long-term temporal trends. Spatial pattern analysis 

were developed to provide an understanding of potential relations between the typical values of water-quality constituents 

and the physical setting of the watershed represented by the 

monitoring station, as well as the relations between observed 

temporal trends and respective watershed characteristics.

Water-Quality Data Structure and Analysis

Water-quality sampling for the SWAS program was 

conducted by program personnel. Grab samples were collected and delivered to the SWAS laboratory at the University 

of Virginia for analysis of chemical constituents. Standard 

laboratory procedures, including the use of stand-alone meters, 

titrations, atomic absorption spectrophotometers, and ion chromatographs, were used as appropriate for the various analyses. 

Details of the method used for each analysis are maintained in 

the NPSTORET database. Sample collection and analysis for 

the SHS study were performed similarly to the SWAS program, although USGS staff were responsible for the collection 

of samples. Water-quality measurements made in association 

with the Shenandoah National Park AM and Fisheries Monitoring Programs did not include any sample collection. For 

these efforts, only field measurements of physical parameters 

(water temperature, pH, specific conductance, and dissolved 

oxygen) were taken. These measurements were made using 

sensors placed in the stream at the time of data collection.

Of the water-quality data compiled in the comprehensive 

database, only the SWAS program data have been collected for 

a sufficient time period and at a sufficient temporal resolution 

to support trend analysis. The dataset used for temporal trend 

analysis of water quality, therefore, consisted of data from 

15 sites (fig. 3) spanning a range of time periods, with the 

longest spanning over 30 years (table 3).

Table 2. Results of principal component analysis on 

environmental factors for monitored watersheds.

[%, percent; <, less than; m, meter; NHD, National Hydrography Dataset]

Percentage of variance 

explained and p-value

Principal 

component 1

Principal 

component 2

Principal 

component 3

Individual 35.544 14.178 12.08

Cumulative 35.544 49.723 61.803

p-value <0.0001 <0.0001 <0.0001

Correlation

Watershed area 0.12 –0.71 0.04

Mean watershed 

elevation

0.85 0.39 0.18

Maximum watershed 

elevation

0.79 –0.20 0.35

Minimum watershed 

elevation

0.41 0.77 –0.02

Precipitation 0.58 –0.10 0.47

Mean watershed slope –0.85 –0.04 0.32

Siliciclastic –0.87 0.38 0.03

Granitic 0.41 –0.27 0.69

Basaltic 0.55 –0.14 –0.74

Carbonate –0.25 0.10 0.03

Developed, open space 0.69 0.25 –0.18

Black locust 0.43 –0.23 –0.68

Open space 0.39 0.40 –0.31

Chestnut oak –0.91 0.20 0.11

Yellow poplar 0.36 –0.67 0.28

Cove hardwood 0.44 –0.39 –0.42

Pine –0.69 0.34 –0.01

Red oak 0.82 0.28 0.01

Hemlock 0.66 0.37 0.33

Aspect –0.26 0.09 –0.34

Channel slope 100 m 0.39 0.23 0.29

NHD reach slope 0.44 0.52 0.02
10 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

N 0 10 20 KILOMETERS

0 5 10 20 MILES

EXPLANATION

!

Shenandoah National Park 

 boundary

Streams

Water-quality trend monitoring 

 station and identifier

Siliciclastic

Bedrock geology

Granitic

Basaltic

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

VT36

PAIN

MAD2

DR01

WOR1

VT53

STAN

VT66

VT75

VT58

NFDR

VT62

VT61

VT51

PINE

Figure 3. Map showing monitoring stations used for temporal trend analysis of water-quality constituents. Location 

of study area provided in figure 1.

Figure 3. Monitoring stations used for temporal trend analysis of water-quality constituents. Location of study 

area provided in figure 1.
Analysis of Water-Quality Data 11

Spatial variability in water-quality conditions throughout the park was evaluated using graphical and statistical 

techniques. Graphical analyses consisted of deriving spatially 

interpolated maps from median constituent values for every 

monitoring station from every monitoring program with such 

data. These analyses provide reasonable information about 

general water-quality patterns throughout the park, because 

sites containing water-quality data are well distributed within 

the park. In general terms, the interpolated maps are intended 

to depict the spatial distribution of water-quality conditions, 

but additional fine-scale variability is probably not represented 

well by these generalizations because these simple interpolations are based on site proximity. Statistical analyses included 

analysis of correlations between water-quality trend slopes and 

median values and the principal component axes for watershed 

characteristics identified earlier herein at those sites where 

trend analyses were performed. These correlation analyses 

were performed using the SWAS program data that were used 

in the analysis of temporal trends.

Spatial patterns in water quality, or associations with 

spatially variable watershed characteristics, were evaluated statistically using correlation analysis (Kendall’s tau) of the waterquality medians and trend slopes with the principal component 

axes that represent the physical characteristics of the monitored 

watersheds. Additionally, scatterplots were generated and bestfit lines were drawn to better characterize the relation between

the principal components and water-quality variables. 

Accounting for Streamflow in Water-Quality 

Trend Assessments 

Variability in streamflow often is chiefly responsible for 

variability in measurements of water quality, as water contributions from various sources may be masked or enhanced 

depending on the flow regime. Changes in streamflow can 

result in the dilution or concentration of water-quality constituents, and removal of this variability can improve the power 

to detect trends (Schertz and others, 1991). Because of these 

dependencies on hydrologic conditions, it was necessary to 

characterize the hydrologic conditions that existed at the time 

of sample collection. Although some of the data-collection 

efforts within the park included measurements of streamflow 

or some general characterization of flows or water levels, the 

water-quality samples used for trend analysis did not. Therefore, it was necessary to develop a representative characterization of flow conditions for all sampling events across the park. 

The SWAS program operates streamgages on five streams 

in the park to generate long-term continuous timeseries 

records of streamflow (table 4). These five streamgages represent a range of watershed sizes and geologic characteristics 

(table 4), and therefore, differing distributions of streamflow 

(fig. 4) believed to collectively represent hydrologic conditions 

throughout the park.

Daily mean streamflow data from the five SWAS 

streamgages were used to compute an approximation of 

hydrologic conditions across the park. The approximation 

was developed by first evaluating the distribution of measured 

daily mean streamflow values available for the period of 

record for each gaged site. From the cumulative distribution 

frequency (fig. 4) for each site, the non-exceedance probability was computed for each daily mean streamflow value. The 

mean of the collection of non-exceedance probabilities for 

each day in the study period was then computed and used as a 

representation of hydrologic conditions, as needed for the statistical analyses of water quality and aquatic biota that follow. 

To confirm the approximation just described was representative of actual hydrologic conditions, overall mean 

non-exceedance values were compared to site-specific 

non-exceedance values and site-exclusive non-exceedance 

values (overall mean non-exceedance computed excluding 

the site of interest). Additionally, concentration-streamflow 

relations were compared with concentration-streamflow relations using the overall mean non-exceedance probabilities, 

the site-specific non-exceedance probabilities, and the siteexclusive non-exceedance probabilities instead of streamflow. 

Table 3. Shenandoah Watershed Study (SWAS) monitoring 

stations, periods of record, years of record, and number of waterquality results contained in the dataset used for trend analysis.

Site 

identifier Site name Period of record

Years of 

record

Number of 

results

DR01 Deep Run 11/02/1979–

12/20/2010

31.2 25,842

MAD2 Madison Run 

(upper)

03/31/1981–

07/26/2010

29.3 10,835

NFDR North Fork Dry 

Run

03/12/1987–

12/21/2010

23.8 23,107

PAIN Paine Run 08/14/1987–

12/20/2010

23.4 30,265

PINE Piney River 10/21/1995–

12/21/2010

15.2 26,172

STAN Staunton River 09/04/1996–

12/21/2010

14.3 18,909

WOR1 White Oak Run 11/02/1979–

12/20/2010

31.2 26,896

VT36 Meadow Run 08/16/1987–

10/28/2010

23.2 2,041

VT51 Jeremy’s Run 04/25/1987–

10/28/2010

23.5 1,923

VT53 Twomile Run 08/14/1987–

10/28/2010

23.2 2,085

VT58 Brokenback 

Run

08/12/1987–

10/29/2010

23.2 2,005

VT61 N.F. Thornton 

River

04/26/1987–

10/28/2010

23.5 1,923

VT62 Hazel River 04/26/1987–

10/28/2010

23.5 1,934

VT66 Rose River 04/25/1987–

10/29/2010

23.5 1,935

VT75 White Oak 

Canyon Run

07/31/1990–

10/29/2010

20.3 1,728
12 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Correlations of streamflow non-exceedance probabilities 

between gaged sites indicated that, although the magnitude of 

flow varied among sites, the non-exceedance probabilities of 

the daily values were strongly correlated, indicating that the 

mean non-exceedance probabilities can be used as a representative approximation of hydrologic conditions throughout the 

park (table 5).

Relations between the computed flow statistics and 

ANC were evaluated using linear regression to evaluate the 

appropriateness of the statistic for flow normalization of 

water-quality data. Regressions of ANC and the overall mean 

non-exceedance probabilities, the site-specific non-exceedance 

probabilities, and the site-exclusive non-exceedance probabilities were evaluated to determine if the use of the mean nonexceedance probability effectively reproduced the site-specific 

concentration-streamflow relation. Consistency across the 

models in the amount of variability in ANC explained by the 

flow statistics, determined from the coefficient-of-determination (R2

) and the root mean squared error (RMSE), demonstrated that the mean probability describes variability in waterquality constituents at a level consistent with the site-specific 

relations (table 6). Furthermore, overlap of the 95-percent 

confidence intervals (CI) for the intercepts and slopes of the 

models for 4 of the 5 sites demonstrates that regression lines 

from each of the flow statistics are generally not significantly 

different (table 6).

Temporal Trends in Water Quality

Temporal trends in water quality were evaluated using 

the Seasonal-Kendall trend test (Hirsch and Slack, 1984) as 

programmed in S-ESTREND (Slack and others, 2003), a 

version of EStimateTREND (ESTREND; Schertz and others, 

1991) written in the S programming language (Venables and 

Ripley, 2000). The Seasonal-Kendall test is a non-parametric 

test of monotonic trend based upon the Mann-Kendall test 

(Mann, 1945; Kendall, 1975) that accounts for seasonal variability by comparing data collected within the same season. 

The Seasonal-Kendall trend test was used because of its ability 

to efficiently handle multiple statistical challenges common to 

trend analysis of environmental data, including non-normality, 

seasonality (serial correlation), streamflow dependence, and 

the violation of statistical assumptions for parametric analysis 

(Schertz and others, 1991; Zipper and others, 2002).

Trend analyses of water-quality constituents pertaining to 

acidic deposition and acidification of stream water, including 

pH, ANC, NO3

–

, and SO4

2–, were performed. Trend analyses 

were conducted at all monitoring stations within the park for 

which the minimum criteria of 5 years of data collection with 

a sample count greater than 10 (Schertz and others, 1991) 

was satisfied. Time periods of 5, 10, 14, 20, and 30 years, all 

ending in 2009, were used to evaluate short-term (recent) to 

long-term trends in the constituents of interest. 

Table 4. Watershed area, geology, and period of record for five streamgages in Shenandoah 

National Park used to determine streamflow statistics.

[km2, square kilometer]

Station 

identifier Stream name

Watershed 

area 

(km2)

Predominant 

geology

Record begin 

date

Record 

end date

WOR1 White Oak Run 5.1 Siliciclastic 11/01/1979 01/17/2011

NFDR North Fork Dry Run 2.3 Granitic 10/01/1987 03/01/2011

STAN Staunton River 10.7 Granitic 09/02/1992 02/14/2011

PINE Piney River 12.4 Basaltic 09/30/1992 02/15/2011

PAIN Paine Run 12.7 Siliciclastic 10/01/1992 02/07/2011

Station identifier

EXPLANATION

NFDR

PAIN

PINE

STAN

WOR1

Non-exceedance probability

0

0.25

0.50

0.75

1.00

0.001 0.01 0.1 1 10 100

Streamflow, in cubic feet per second

Figure 4. Graph showing cumulative distribution frequency 

of daily mean streamflow from five stream gages used in 

Shenandoah Watershed Study in Shenandoah National Park.

Figure 4. Cumulative distribution frequency of daily mean 

streamflow from five streamgages used in Shenandoah 

Watershed Study in Shenandoah National Park.
Analysis of Water-Quality Data 13

Table 5. Correlation coefficients (Pearson’s r) for non-exceedance probabilities from the five 

Shenandoah Watershed Study (SWAS) streamgages and the mean non-exceedance probability 

computed from those five streamgages. 

[All coefficients significant for α=0.01. —, value not shown; N/A, not applicable]

Station 

identifier or 

statistic

NFDR PAIN PINE STAN WOR1 Mean

Siteexclusive 

mean

NFDR 1.000 — — — — — 0.870

PAIN 0.857 1.000 — — — — 0.935

PINE 0.880 0.888 1.000 — — — 0.923

STAN 0.832 0.862 0.868 1.000 — — 0.893

WOR1 0.831 0.927 0.860 0.844 1.000 — 0.885

Mean 0.937 0.959 0.952 0.932 0.962 1.000 N/A

Table 6. Statistics for regression models using various flow statistics to estimate acidneutralizing capacity, by site.

[ANC, acid-neutralizing capacity; RMSE, root mean squared error; R2, coefficient of determination; %, percent; CI, 

confidence interval]

Non-exceedance 

probability

Mean 

ANC RMSE R2

Slope Intercept

Lower 95% 

CI

Upper 95% 

CI

Lower 95% 

CI

Upper 95% 

CI

NFDR

Site exclusive 63.48 25.58 0.28 –62.85 –52.17 89.73 95.92

Mean 63.48 25.30 0.30 –65.34 –54.63 90.75 96.91

Site specific 63.34 26.21 0.24 –55.71 –45.08 85.43 91.55

PAIN

Site exclusive 6.15 3.92 0.56 –17.02 –15.20 13.82 14.87

Mean 6.15 3.86 0.57 –17.17 –15.37 13.93 14.97

Site specific 6.15 3.87 0.57 –16.41 –14.67 13.61 14.63

PINE

Site exclusive 233.05 41.68 0.67 –230.87 –210.72 338.90 350.46

Mean 233.05 41.19 0.68 –231.09 –211.25 338.89 350.26

Site specific 233.15 44.16 0.63 –208.45 –188.65 326.39 337.84

STAN

Site exclusive 84.50 15.00 0.35 –43.37 –36.43 102.76 106.77

Mean 84.50 14.97 0.35 –43.80 –36.83 102.99 107.01

Site specific 84.40 15.16 0.31 –38.48 –31.74 100.33 104.29

WOR1

Site exclusive 28.37 14.63 0.44 –53.32 –46.43 53.25 57.37

Mean 28.89 13.47 0.50 –56.53 –50.82 56.64 60.09

Site specific 28.89 13.46 0.50 –54.91 –49.37 56.06 59.44
14 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Some of the variability in water quality often can be 

attributed to seasonal variations resulting from phenomena 

such as biological activity, changes in constituent source, or 

changes in the predominant water source, and such variability 

can mask temporal trends if not removed (Hirsch and others, 1982; Schertz and others, 1991). The Seasonal-Kendall 

trend test reduces the seasonal effect by grouping the data into 

seasons and then comparing data collected within the same 

season (Hirsch and Slack, 1984). The selection of the number 

of seasons to use in a given trend test was based on either 

(1) an understanding of seasonal variability in the constituent of interest, or (2) the number of annual samples collected 

according to a consistent sampling frequency (Schertz and others, 1991). For this study, the sampling frequency was generally biweekly, monthly, or quarterly. Because of the potentially 

confounding effects of serial correlation among frequently 

collected samples, and because seasonal differences rarely 

occur on a scale finer than 1 month, the S-ESTREND program 

does not allow for seasonal definitions finer than 12 per year 

(monthly; Schertz and others, 1991). Consequently, quarterly seasonal definitions were used for all sites, and, where 

possible, monthly seasonal definitions were used in order to 

compare the effects of different sampling frequencies on the 

ability to detect trends.

Trend analyses were conducted on both constituent 

concentrations and streamflow-adjusted concentrations. 

Streamflow adjustment was accomplished using the stagewise 

method (Hirsch and others, 1982) in which a regression of 

streamflow versus concentration is developed and the trend 

test is performed on the residuals of the regression. In this 

case, the residual equals the observed concentration minus 

the predicted concentration. The flow-concentration regression was performed in this study using LOESS models, a 

procedure similar to LOcally WEighted Scatterplot Smoothing (LOWESS). LOESS is a robust approach for flow adjustment when regression assumptions are violated and (or) when 

thorough evaluation of model fit is not feasible because of 

the large number of possible models generated when working with numerous constituent/site combinations (Schertz 

and others, 1991; Slack and others, 2003). Streamflow data 

were not available for all monitoring stations where waterquality data were collected; therefore, estimations of hydrologic conditions, as described in the previous section, were 

used in the flow-adjusted trend models. All default settings 

for seasonal periods, minimum data requirements, and flow 

adjustment within the S-ESTREND program were utilized. 

Alpha (α)=0.10 was specified as the cutoff value for statistical 

significance of the trend slopes. 

Trend analyses were furthered by evaluating the significance of trends park-wide and grouped by major geologic 

classification. The intent of this analysis was to evaluate the 

potential for geologic effects on water-quality trends and 

to provide a general overview of the direction of change in 

water-quality constituents throughout the park. This was 

accomplished using a 2-sided one-sample t-test on all trend 

slopes for a given constituent and time period. The t-test is 

used to test the null hypothesis that the trend slope is equal 

to zero. A significant result from the t-test indicates that the 

mean trend slope is significantly different from zero, thereby 

indicating the presence of a geology class-wide trend, regardless of the significance of the individual trend results. Similar 

analyses to evaluate the potential effects of watershed area, 

given the importance of this factor in the PCA, were intended 

but not completed because of disparities in the distributions of 

watershed areas within geologic classes.

Determination of Appropriate Trend-Model 

Formulations

Comparisons of slopes for concentration trends with 

flow-adjusted concentration trends indicated that the two 

approaches yielded consistent patterns for ANC, pH, and 

SO4

2– concentrations. That is, the magnitude and significance 

of trend slopes were similar whether concentration or flowadjusted concentration was used in the analyses (fig. 5). Therefore, subsequent analyses focus on concentration-only trends 

for ease of interpretation and because trends in concentration 

are more pertinent in evaluating the effects of water quality on 

aquatic biota.

The SWAS program has collected water-quality data at 

various frequencies across sites and over time. As described 

earlier, only the monthly and quarterly data are appropriate for 

use in analyses of temporal trends. To evaluate the effect of 

sampling frequency on the results of temporal trend analysis, 

the trend analyses were computed at the sites where monthly 

data were collected using both monthly and quarterly seasonal 

definitions. This evaluation was accomplished by comparing the statistical significance of the computed concentration 

trend slopes using the 4-season and 12-season definitions 

(fig. 6). Generally, the significance (α=0.1) of the computed 

trend slope was similar between the 4-season and 12-season 

definitions, and the use of the 4-season definition allows trend 

analyses to be conducted at a greater number of sites; therefore, the temporal trends in concentration (not flow adjusted) 

presented hereafter are based upon the 4-season definition 

(quarterly sampling).
Analysis of Water-Quality Data 15

Temporal Trends in Water-Quality 

Constituents

Trends in ANC, pH, and SO4

2– were analyzed 

using methods just detailed. Evaluation of trends 

in NO3

–

 was unsuccessful because of the high proportion of censored results; greater than 40 percent 

of NO3

–

 results were below the reporting limit of 

the laboratory analysis. The statistical methods 

employed within S-ESTREND can effectively 

detect trends in datasets having a small proportion of censored values, but once that proportion 

exceeds about 5 percent, the results are impaired 

(Schertz and others, 1991).

Temporal trends in ANC were strongly correlated with median ANC values, which were 

strongly related to the underlying geology of the 

watershed (fig. 7). These relations were most 

obvious when comparing trend slopes and median 

ANC values between siliciclastic and basaltic 

watersheds. Most negative trend slopes in ANC 

were observed in siliciclastic watersheds, which 

also had the lowest median ANC values. These 

patterns suggest that streams draining siliciclastic 

watersheds were not only more acidic than streams 

in watersheds of other geologic types, but also 

showed stronger evidence of being increasingly 

affected by acidic deposition over time, especially 

in the more recent time periods (5 and 10 years). 

Although these negative trend slopes are small, 

on the order of –1 μeq/L per year, these results 

may be biologically relevant because median ANC 

values in these streams are only about 12 μeq/L, 

making the proportional change per year relatively 

large. Furthermore, decreasing ANC trends were 

statistically significant geology-class-wide only 

for streams in siliciclastic watersheds (for the 5-, 

10-, and 14-year periods), and class-wide significant increasing ANC trends were observed only in 

basaltic watersheds (for the 20-year period). Longterm significant trends and short-term non-significant trends in the basaltic watersheds indicate 

improving conditions, because ANC is increasing 

in these watersheds. The degree of improvement 

is relatively small, however, because the median 

representative value of ANC in these watersheds 

is nearly 200 μeq/L and the significant long-term 

trends have slopes of 1 to 2 μeq/L per year. ANC 

trends in the granitic watersheds show a distinct 

pattern of change, with trends split proportionately between increasing and decreasing slopes, 

only two of which are statistically significant. No 

class-wide significant ANC trends were observed 

in granitic watersheds. 

EXPLANATION

Trend slope

Significant (α = 0.1) 

 flow adjusted

Not-significant (α = 0.1) 

 flow adjusted

Significant (α = 0.1) 

 concentration

Not-significant (α = 0.1) 

 concentration

Figure 5. Graphs showing relations between 

concentration trend slope and flow-adjusted 

trend slope for A, acid neutralizing capacity 

(ANC), B, pH, and C, sulfate (SO4

2-).

–5

–4

–3

–2

–1

0

1

2

3

4

5

–5 –4 –3 –2 –1 0 1 2 3 4 5

–0.05

–0.04

–0.03

–0.02

–0.01

0

0.01

0.02

0.03

0.04

0.05

Flow-adjusted trend slope

–0.05 –0.04 –0.03 –0.02 –0.01 0 0.01 0.02 0.03 0.04 0.05

–3.0

–2.5

–2.0

–1.5

–1.0

–0.5

0

0.5

1.0

1.5

2.0

2.5

3.0

A

B

C

–3.0 –2.5 –2.0 –1.5 –1.0 –0.5 0 0.5 1.0 1.5 2.0 2.5 3.0

Concentration trend slope

Figure 5. Relations between concentration trend slope and 

flow-adjusted trend slope for A, acid neutralizing capacity (ANC), 

B, pH, and C, sulfate (SO4

2–).
16 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Trends in pH were more variable and less indicative of 

relations with underlying geology than ANC trends (fig. 8). 

Long-term trends (10-, 14-, and 20-years) across all geologic 

categories generally indicated increasing pH values, or recovery from acidification, with significant class-wide trends in the

granitic watersheds (10- and 20-year trends), basaltic watersheds (10-, 14-, and 20-year trends), and siliciclastic watersheds 

(14-year trends). Significant long-term trends were detected for

trend slopes of less than 0.01 to more than 0.03 unit. Although 

these trend slopes are small, they equate to increases of 0.2 

to 0.3 pH unit over the 10- to 30-year duration of the periods 

evaluated, which could be biologically meaningful. 

Short-term trends (5 year) in pH exhibit uncharacteristic 

declines compared to the longer-term trends, and these trends 

were significant class-wide for each of the three geologic 

classes. Short-term trends are highly susceptible to the effects of 

climatic variability and were probably affected by dry conditions 

early in the trend period, followed by wet conditions later in 

the period. Such climatic conditions would result in a decreasing trend in pH as greater volumes of acidic deposition were 

delivered to the watershed during the latter part of the 5-year 

period. This effect also is apparent in the short-term ANC trends 

in siliciclastic watersheds, where ANC values are very low, but 

was not detected in the basaltic and granitic watersheds where 

background ANC is much higher, exemplifying the sensitivity of 

the streams underlain by siliciclastic geology. 

Trends in SO4

2– show variability similar to the trends in 

pH (fig. 9). Significant long-term SO4

2– trends in the siliciclastic and granitic watersheds had slopes of approximately 

0.5 μeq/L per year, which yields a long-term increase of 

approximately 10 μeq/L over the 20-year period. Only the 

14- and 20-year trends in the granitic class, however, yielded 

a significant increasing class-wide trend. In contrast, longterm SO4

2– trends in the basaltic watersheds revealed declining 

SO4

2– concentrations in those watersheds, with trend slopes as 

great as approximately –1 μeq/L per year, yielding an approximately 20-μeq/L decline over the 20-year period. The 10-year 

trend results for the basaltic class yielded a significant classwide decreasing trend, as did the 14-year trend in the siliciclastic class. Although the short-term trends in SO4

2– exhibit 

greater variability than the long-term trends, the overall pattern indicates increasing SO4

2– concentrations in the granitic 

watersheds, decreasing SO4

2– concentrations in the basaltic 

watersheds, and mixed results in the siliciclastic watersheds. 

With the exception of the highly variable 5-year trends, these 

results appear to indicate decreasing trends in recent years.

Acid neutralizing capacity (ANC) pH Sulfate (SO4

2-)

p-value of 12-season trend

0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0

p-value of 4-season trend

EXPLANATION

5 years

10 years

20 years

30 years

1:1 line

Period of trend 

 analysis

p = 0.1

Figure 6. Graphs showing trend-slope p-values for 4-season sampling definition against p-values for 12-season 

sampling definition.

Figure 6. Trend-slope p-values for 4-season sampling definition against p-values for 12-season sampling definition.
Analysis of Water-Quality Data 17

Figure 7. Graphs showing acid neutralizing capacity (ANC) trend slopes and significance by geologic 

category and analysis period. Years with significant (α = 0.10) class-wide trends are denoted 

with an asterisk.

300

250

200

150

100

50

0

Median ANC, in 

 microequivalents 

 per liter

EXPLANATION

Not significant 

 (α = 0.10)

Trend slope

Significant 

 (α = 0.10)

Granitic Basaltic Siliciclastic

Trend slope, in microequivalents per liter per year

–5

–4

–3

–2

–1

0

1

2

3

4

5

5 10 14 20 30 5 10 14 20* 30 5* 10* 14* 20 30

Period of analysis, in years

Figure 7. Acid neutralizing capacity (ANC) trend slopes and significance by geologic category and 

analysis period. Years with significant (alpha = 0.10) class-wide trends are denoted with an asterisk.

7.50

7.00

6.50

6.00

5.50

Median pH

EXPLANATION

Trend slope

Not significant 

 (α = 0.10)

Significant 

 (α = 0.10)

Granitic Basaltic Siliciclastic

–0.05

–0.04

–0.03

–0.02

–0.01

0

0.01

0.02

0.03

0.04

5* 10* 14 20* 30 5* 10* 14* 20* 30 5* 10 14* 20 30

Trend slope, in microequivalents per liter per year

Period of analysis, in years

Figure 8. pH trend slopes and significance by geologic category and analysis period. Years with 

significant (alpha = 0.10) class-wide trends are denoted with an asterisk.
18 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Practical Significance of Trend Analyses

It is often desirable to statistically evaluate the power 

of a given analysis, or to determine what the smallest detectable change is for a given sampling design. Although a formal 

power analysis was not performed, insight may be gained as to 

the power of these trend analyses by evaluating the practical 

significance of the smallest detected significant trends for each 

constituent.

Significant trends in ANC were detected in 15 of the 

62 trend analyses performed (fig. 7). The mean absolute value 

of the significant trend slopes for ANC was approximately 

0.9 μeq/L per year, whereas the minimum absolute value 

of the significant trend slopes was approximately 0.2 μeq/L 

per year. These trends slopes were smallest at sites having 

median long-term ANC values less than 5 μeq/L, or those sites 

most affected by acidification. These results indicate that the 

quarterly sampling scheme produces sufficient data to detect 

long-term ANC trends of practical significance. 

Significant trends in pH were detected in 17 of the 

62 trend analyses performed (fig. 8). The mean absolute value 

of the significant trend slopes for pH was approximately 

0.02 units per year, whereas the minimum absolute value of 

the significant trend slopes was approximately 0.003 units per 

year. The trend slopes were smallest at sites having the lowest 

median long-term pH, or those sites that are most degraded in 

terms of pH. The magnitude of these trends is extremely small 

and probably well below the level of practical significance. It 

is valuable to be able to detect such small changes over time, 

however, because minor changes in pH may ultimately lead to 

the crossing of toxic thresholds for aquatic biota, in either the 

recovery or degradation direction (for example, see Deviney 

and others, 2012). These results indicate that the quarterly 

sampling scheme produces sufficient data to detect long-term 

pH trends of practical significance, particularly at the most 

degraded sites. 

Significant trends in SO4

2– were detected in 7 of the 

62 trend analyses performed (fig. 9). The mean absolute value 

of the significant trend slopes for SO4

2– was approximately 

0.6 μeq/L per year, whereas the minimum absolute value 

of the significant trend slopes was approximately 0.2 μeq/L 

per year. These trend slopes all represent trends of less than 

1 percent per year when evaluated in terms of the median 

SO4

2– concentration. The detection of such small changes in 

SO4

2– indicate that the quarterly sampling scheme produces 

sufficient data to detect long-term SO4

2– trends of practical 

significance. 

EXPLANATION Granitic Basaltic Siliciclastic

–3

–2

–1

0

1

2

3

5 10 14* 20* 30 5 10* 14 20 30 5 10 14* 20 30

Trend slope

Not significant 

 (α = 0.10)

Significant 

 (α = 0.10)

30

40

50

60

70

80

90

100

110

120

Median sulfate (SO4

2–

) 

 concentration, 

 in microequivalents 

 per liter

Figure 9. Graphs showing sulfate (SO4

2-) trend slopes and significance by geologic category and 

analysis period. Years with significant (α = 0.10) class-wide trends are denoted with an asterisk.

Trend slope, in microequivalents per liter per year

Period of analysis, in years

Figure 9. Sulfate (SO4

2–) trend slopes and significance by geologic category and analysis period. 

Years with significant (alpha = 0.10) class-wide trends are denoted with an asterisk.
Analysis of Water-Quality Data 19

Temporal Trends in Water Temperature

Water-temperature trends, particularly warming trends, 

in the park are of interest because of the potential implications for cool- and cold-water organisms, such as brook trout, 

that are native to park waters and prefer water temperatures 

below 17 °C (Cherry and others, 1977). Significant increases 

in stream-water temperature were detected at numerous sites 

in the park, predominantly during the 10-year period (fig. 10). 

Increasing trend slopes in the annual median, mean, and maximum water temperature were often large, ranging from 0.14 to 

0.76 °C per year, with a median of 0.32 °C per year. The longer-term trend results (14-year and 30-year) also indicate rising water temperatures, particularly annual maximums. Only 

one site (WOR1) had sufficient data to support evaluation of 

the 30-year trend in water temperature, and the results for that 

site indicated a small but statistically significant increase in 

annual mean, median, and maximum water temperature for 

the 30-year period ending in 2009. The trend slope for mean 

annual water temperature over this period was approximately 

0.04 °C per year, which equates to an increase of 1.2 °C over 

the 30-year period. As with the 5-year trends in water quality, 

the 5-year trends in water temperature were highly variable.

One explanation for the observed increases in streamwater temperatures in the park is that they were due to 

increased solar radiance associated with forest defoliation 

from gypsy moth infestation. For this hypothesis to be consistent with the data, water temperatures would be expected to 

increase during the period of defoliation, followed by decreasing water temperatures, as canopy cover was reestablished and 

streams were once again shaded. In many cases, increasing 

water temperatures were observed prior to, during, and after 

the period of most intense defoliation in the park (1990–95; 

fig. 11). Thus, the observed patterns do not support the hypothesis that gypsy moth defoliation is responsible for increased 

water temperatures.

An alternative hypothesis for observed increases in 

water temperatures in the park is that air temperature has 

increased in association with climate change. In order to test 

the climate-change hypothesis, changes in water temperature 

were evaluated over larger spatial and longer temporal scales. 

For these analyses, long-term water-temperature records were 

obtained from nine USGS streamgages located within and near 

the park (fig. 12). These data represent a range of watershed 

areas, from small watersheds similar in size to those in the 

park, to much larger watersheds that incorporate the drainages 

of the park. The period of record for these sites is sufficiently 

long to examine trends at a temporal scale relevant to climate 

change (table 7). White Oak Run, one of the streams for which 

water-temperature trends were computed from the SWAS 

dataset, also had a USGS streamgage that was operated at a 

separate location. Measurements at the White Oak Run USGS 

streamgage were made independently from the SWAS efforts, 

and therefore represent an independent verification of the 

water-temperature trends observed in the SWAS dataset.

EXPLANATION

Not significant (α = 0.10)

Significant (α = 0.10)

Trend slope, in degrees Celsius per year

–0.5

0

0.5

1.0

–0.5

0

0.5

1.0

–0.5

0

0.5

1.0

1.5

1.5

1.5

1.5

–0.5

–1.0

–1.0

–1.0

–1.0

0

0.5

1.0

Maximum

Mean

Median

Minimum

5 10 14 20 30

Period of analysis, in years

Figure 10. Graphs showing trend slopes and significance 

of annual water-temperature characteristics by analysis 

period.

Figure 10. Trend slopes and significance of annual 

water-temperature characteristics by analysis period.
20 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

EXPLANATION

Mean

Minimum

Maximum

Water temperature, in degrees Celsius

0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

0

5

10

15

20

25

30

DR01 NFDR PAIN

PINE STAN VT36

VT51 VT53 VT58

VT61 VT62 VT66

VT75 WOR1

1980 1985 1990 1995 2000 2005 2010 1980 1985 1990 1995 2000 2005 2010 1980 1985 1990 1995 2000 2005 2010

Year

Figure 11. Graphs showing time series of mean, minimum, and maximum annual water temperature for 14 sites, with period of 

defoliation indicated by vertical dotted lines.

Figure 11. Time series of mean, minimum, and maximum annual water temperature for 14 sites, with period of defoliation 

indicated by vertical dotted lines.

Table 7. Watershed area, mean watershed elevation, and period of record for U.S. Geological Survey (USGS) 

streamgages used for water-temperature trend analysis.

[km2, square kilometer; m, meter]

Station 

identifier Station name

Watershed 

area 

(km2)

Mean watershed 

elevation 

(m)

Period of 

record

Years of 

record

01620500 North River Near Stokesville, Va. 45 972 1941–2011 70

01621050 Muddy Creek At Mount Clinton, Va. 37 617 1993–2011 18

01627500 South River At Harriston, Va. 548 756 1925–2011 86

01628060 White Oak Run Near Grottoes, Va. 5 709 1978–1996 18

01628500 S.F. Shenandoah River Near Lynnwood, Va. 2,795 826 1930–2011 81

01629500 S.F. Shenandoah River Near Luray, Va. 3,554 781 1925–2011 86

01631000 S.F. Shenandoah River At Front Royal, Va. 4,233 743 1930–2011 81

01662800 Battle Run Near Laurel Mills, Va. 67 466 1958–2011 53

01665500 Rapidan River Near Ruckersville, Va. 297 650 1989–2011 22
Analysis of Water-Quality Data 21

Air-temperature records from two meteorological stations 

within the park (Luray 5E and Big Meadows) were used to 

evaluate long-term trends in air temperature. These air-temperature data were retrieved from the National Climatic Data 

Center website (National Oceanic and Atmospheric Administration, 2012) as monthly mean air temperature, monthly mean 

maximum air temperature, and monthly mean minimum air 

temperature. The monthly mean values were averaged over 

each calendar year to generate annual mean air temperatures, 

annual mean minimum air temperatures, and annual mean 

maximum air temperatures for use in the trend analyses. Data 

included records for the period 1941–2011 at the Luray site 

and 1935–2011 at the Big Meadows site; however, recent 

data were incomplete for both sites, so the trend analysis for 

each was performed only through 2001. Trend analyses for 

air temperature were performed using simple linear regression (1) for the entire period of the analysis dataset to evaluate 

long-term trends, and (2) from 1979 to the end of the period of 

analysis dataset for comparison with the maximum period of 

record for water temperature in SHEN streams.

Significant water-temperature trends of similar magnitude to those observed in the park were detected at six of the 

nine USGS streamgages in and around the park (fig. 13). Significant water-temperature trends were most common at stations having a period of record similar to the period of record 

of the SWAS sites, and at stations with the smallest watershed 

area (table 7; fig. 13). Furthermore, water-temperature trends 

observed at the USGS streamgage on White Oak Run near 

Grottoes, Va., over a period of record that overlaps the SWAS 

dataset are of similar magnitude to the water-temperature 

!

!

#

#

#

#

#

#

#

#

#

Luray 5 E

Big Meadows

01665500

EXPLANATION

Shenandoah National Park

# USGS streamgages and 

 identifier

! NOAA weather station and 

 identifier

Luray 5 E

01665500

01662800

01631000

01629500

01628500

01628060

01627500

01621050

01620500 River

Rapidan

Robinson

River

Hazel

Thornton

River

Goose

Creek

Cedar

Creek

North

Fork

Shenandoah

Shenandoah

River

Shenandoah

Fork

South

River

River

Dry

Middle

River

River

South

79° 78°

39°

38°

Base from U.S. Geological Survey digital 

line graph, 1:2,000,000, 1987

Figure 12. Map showing U.S. Geological Survey (USGS) streamgages used in water-temperature trend analysis. 

Location of study area provided in figure 1. NOAA is National Oceanic and Atmospheric Administration.

0 10 20 Miles

0 10 20 Kilometers

VIRGINIA

WEST VIRGINIA

Figure 12. U.S. Geological Survey (USGS) streamgages used in water-temperature trend analysis. Location of study area 

provided in figure 1. NOAA is National Oceanic and Atmospheric Administration.
22 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

trends observed in the SWAS dataset (fig. 13). Lastly, a more 

moderate increasing trend was observed over the much longer 

period of record at the larger South Fork Shenandoah River 

near Lynnwood, Va. streamgage compared to the trends 

observed in the park, suggesting that increases in water temperature were not limited to the small mountainous watersheds 

such as those in the park. The increasing water-temperature 

trends observed in the small, mountainous, streams of the park 

(fig. 10) exceed the magnitude of water-temperature trends 

documented by other studies using similar methods (Kaushal 

and others, 2010). The lesser magnitude trend at the larger 

Lynnwood, Va. Streamgage, compared to the trends observed 

in the park streams, is consistent with those reported in the 

literature, which tend to focus on larger river systems (Kaushal 

and others, 2010).

Trend slope

EXPLANATION

Not significant

Significant

Trend slope

–0.2

–0.3

–0.3

–0.3

–0.1

0

0.1

0.2

0.3

0.4

–0.2

–0.1

0

0.1

0.2

0.3

0.4

–0.2

–0.1

0

0.1

0.2

0.3

0.4

White Oak Run near Grottoes, VA Muddy Creek at Mount Clinton, VA North River near Stokesville, VA

Battle Run near Laurel Mills, VA Rapidan River near Ruckersville, VA South River at Harriston, VA

S.F. Shenandoah River near Lynnwood, VA S.F. Shenandoah River near Luray, VA S.F. Shenandoah River at Front Royal, VA

Maximum

Mean

Median

Minimum

Maximum

Mean

Median

Minimum

Maximum

Mean

Median

Minimum

Water temperature statistic

Figure 13. Annual water-temperature trend slopes and significance for nine U.S. Geological Survey (USGS) streamgages. 
Analysis of Water-Quality Data 23

Comparison of Water-Temperature Trends and 

Air-Temperature Trends

 The trend analysis of air-temperature records indicates that 

air temperature has generally declined over the entire period of 

record (1935–2001; fig. 14). Air-temperature data from 1979 to 

2001, the period for which water-temperature trends were evaluated in the park, however, show increasing trends comparable to 

water-temperature trends observed at sites throughout the park. 

These patterns suggest that increases in water temperature were 

associated with increases in air temperature. 

The longer period of record for the air-temperature data 

provides an important perspective for the water-temperature 

trends observed in the park during recent decades. The period 

of warm air temperatures observed early in the air-temperature 

records (until about 1960; fig. 14) is comparable to the warm 

air temperatures observed in recent decades (1990s–2000s), 

and this pattern also is evident in the long-term water-temperature record at the South Fork Shenandoah River near 

Lynnwood, Va. streamgage. These patterns do not appear to 

be consistent with the climate-change hypothesis because 

increases of similar magnitude occurred much earlier, before 

climate change effects would be expected. It is possible that 

the increasing air temperature trends observed during the earlier decades (before the 1950s), however, were due to timber 

harvesting activities that were widespread during the early part 

of the century and usually involved complete removal of the 

forest (Clarkson, 1964; DeVivo, 1986); and that the subsequent cooling trends observed between 1960 and 1980 were 

associated with the forest re-growth that followed. In fact, 

as late as 1941, more than 15 percent of the park remained 

deforested with nearly 50 percent of the remaining forest 

comprising younger trees (1- to 20-year age class) (Berg and 

Moore, 1941). Therefore, much of the park landscape was still 

recovering from timber harvest and agriculture well into the 

middle of the century. Thus, climate change remains a plausible hypothesis to explain more recent (1980–2010) increasing temperature trends. 

Air-temperature patterns during the last 30 years, a 

period when forest growth has far exceeded forest removal 

throughout the eastern United States (Foster and others, 2004), 

showed increases as large as those observed during the early 

part of the century when much of the Appalachian forest was 

harvested. Other factors, such as the North Atlantic oscillation 

(Visbeck and others, 2001), and increasing urbanization (Kalnay and Cai, 2003), however, also may have contributed to the 

increases in air temperature observed in recent decades in and 

around the park. More research is needed, and it may not be 

possible to determine the precise effect of climate change on 

air temperatures in the park for any particular period of time. 

Temperature trends in air and surface waters observed in this 

study, and similar patterns observed by others (for example, 

Foster and others, 2004), however, are of significant interest to 

resource managers. 

Big Meadows Luray 5 E Lynnwood

Temperature, in degrees Celsius

0

5

10

15

20

25

193019401950196019701980199020002010193019401950196019701980199020002010194019501960197019801990200020102013

Year

Mean minimum air temperature

Mean maximum air temperature

Mean air temperature

Mean water temperature

Trend line not significant

Trend line significant

EXPLANATION

Figure 14. Annual air-temperature characteristics at the Big Meadows and Luray meteorological stations, and annual 

mean water temperature at the Shenandoah River at Lynnwood, Virginia, streamgage, with trend lines.
24 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Spatial Patterns in Water Quality

Spatial patterns in water quality throughout the park, 

particularly those associated with acidification, largely follow 

a north-south gradient coincident with geologic variation 

(fig. 15). Sites with pH data are well distributed throughout 

the park (fig. 15A), and the association of lower pH values in 

the areas underlain by siliciclastic geology (fig. 15) and higher 

pH in the areas underlain by granitic and basaltic geology 

is evident (fig. 15A,E). As would be expected, spatial patterns in ANC (fig. 15B) are similar to the patterns observed 

in pH, although the distribution of sites with ANC data is less 

extensive (fig. 15B). SO4

2– data have been collected at fewer 

sites in the park than pH and ANC, making spatial interpolation of SO4

2– patterns less informative (fig. 15C). Although the 

graphical representation of SO4

2– is coarse, it reveals the same 

general patterns of acidity evident for pH and ANC, namely, a 

generally increasing gradient in SO4

2– and decreasing gradients in pH and ANC from north to south, with the siliciclastic 

region having much higher stream SO4

2– concentrations and 

lower pH and ANC than streams in the remainder of the park. 

Specific conductance, a measure of electrical conductivity 

indicative of total ion concentration in water (Hem, 1985), 

has been measured extensively in the park (fig. 15D). Specific 

conductance of park streams is low, as would be expected for 

relatively undisturbed forested headwater streams. Variability 

within this low range of specific conductance values is associated with the relative solubility of the underlying geology, 

with the siliciclastic and granitic watersheds having lower 

conductivity than the watersheds draining basaltic areas.

Each of the median water-quality values (ANC, pH, 

SO4

2–) were significantly correlated with the first principal 

component, ANC and pH having positive correlations and 

SO4

2– having negative correlation. Additionally, the 10- and 

14-year trend slopes for pH and the 10- and 20-year trend 

slopes for ANC were positively correlated with the first 

principal component. Increases in the first principal component represent increases in elevation and decreases in the mean 

watershed slope and the proportion of siliciclastic geology in 

the watershed. Therefore, these results indicate that pH and 

ANC are greatest, and SO4

2– is lowest, in high-elevation watersheds having flat terrain and a small proportion of siliciclastic 

geology. Furthermore, these results indicate that pH and ANC 

have increased over time in such watersheds—one indication 

of recovery from acidification.

A significant negative correlation with principal component two was detected for median ANC, median pH, ANC 

trends for the 10- and 20-year periods, and the pH trend for 

the 10-year period (fig. 16). Increases in principal component 

two are most strongly influenced by decreasing watershed 

area and increasing minimum watershed elevation (specifically, monitoring station elevation). The correlations indicate 

that ANC and pH values are lower, and trends in pH and ANC 

have lower slopes (some negative), in such watersheds.

The only significant correlations with the third principal component were for measures of SO4

2–, specifically the 

median and 10- and 20-year trends. This component largely 

differentiates granitic geology from the other geologic types, 

with increasing component values reflecting greater proportions of granitic bedrock. The correlations between principal 

component 3 and SO4

2– are positive, indicating a tendency for 

increased SO4

2– concentrations in granitic watersheds.

Overall, this analysis demonstrates that water quality in 

watersheds with higher mean elevations, lower land-surface 

gradients, larger watershed areas, and lesser proportions of 

siliciclastic geology are least affected by acidification and tend 

to be improving over time. In contrast, watersheds with greater 

proportions of siliciclastic and granitic geology, with smaller 

watershed areas, and higher minimum watershed elevation 

tend to be more affected by acidification and are experiencing 

continued degradation.
!

!

!

!

!!!

!!

!

!

!!

!

!

!

!

!

!

!

!

!

!

!!

!

!!

!

!

!

!

!

!

!

!

!

!!

!!

!

!

!!

!

!

!!

!

!

!!

!

!

!

!

!!

!

!

!

!

!

!

!

!

!!

!!

!

!

!!

!!

!

!

!!

!

!

!

!

!

!

! !

!

!!

!!

!!!

!

!

!

! !!

!!

!

!

!

!

!

!

!

!!!

!

!

!

!

!

!!

!

!!

!!

!

!

!

!

!

!!!

!

!

!!!

!

! !

!!

!

!

!

!

!!!!!!!!

!!

!!

!

!

!

!

!

!

!

!

!

!

!

! !

! !!! !!!!!

!!!!

!

!

!

!!

!!!!

!!

!

!

!!!

!!

!!

!

!!!!!!

!!!!

!!!!!! !!

!

!

!

! !

!

!

!!

!

!

!!!

!! !!!

!! !!

!

!

!!!!!!

!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!

!

!!!!!!!!

!!

! ! !

!! !!!!!

!

!

!!!

!

!!

!

!

!!

!!

!!

!

!!!!!

!!!!

!!!!!!

!

!!

!

!

!!

!! !!!

!! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!

!

!

!

!!! !!!!!! !! !!!! ! !!!!!!!!!!!!!!

!! !! !!!!!!

!

!

!

!!! !!! !

!! !! ! !! ! ! ! !! !!! ! ! !!!

! ! !

! !

! !

! !!!!! ! !! !!!! ! !! ! ! !! ! ! ! !! ! !! !! !! !! ! ! !!

! !

!

!!!!!!!!!

!!

!!

! !

!

!

!

!

!

!

!

!

! !!!

! !!

! !

!! !!! ! ! !!

!! ! !!

! ! !!

!!! ! !!

!! ! !!

!

EXPLANATION

Shenandoah National  Park boundary

EXPLANATION

Median specific  conductance, in  microsiemens per  centimeterSpecific conductance  sampling site

17 to 23>23 to 28>28 to 34>34 to 39>39 to 45

EXPLANATION

pH sampling siteMedian pH

4.8 to 5.05.1 to 5.55.6 to 6.06.1 to 6.56.6 to 7.07.1 to 7.57.6 to 8.0

EXPLANATION

Median ANC, in  microequivalents  per literANC sampling site

–15 to 20>20 to 47>47 to 77>77 to 103>103 to 129>129 to 155>155 to 174>174 to 201>201 to 240>240 to 278

EXP

SO4

Me

 c in  p

30 >>>>>>>>

Figure 15. Maps showing spatial patterns in 

A, pH, 

B, acid neutralizing capacity (ANC), 

C, sulfate (SO42-), and 

D, specific conductanc

geology, in Shenandoah National Park. Location of study area provided in figure 1.

A

B

C

D E

N

0 10 20 KILOMETERS0 10 20 MILES

Siliciclastic

Bedrock geology

GraniticBasaltic

Figure 15. Spatial patterns in 

A, pH, 

B, acid neutralizing capacity (ANC), 

C, sulfate (SO42–), and 

D, specific conductance, and 

E, underlying 

geology, in Shenandoah National Park. Location of study area provided in figure 1.
26 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

0

100

200

300

Median of acid 

neutralizing capacity

5-year trend of acid 

neutralizing capacity

10-year trend of acid 

neutralizing capacity

14-year trend of acid 

neutralizing capacity

20-year trend of acid 

neutralizing capacity

Median of pH

5-year trend of pH

10-year trend of pH

14-year trend of pH

20-year trend of pH

Median of sulfate 

concentration, in 

microequivalents 

per liter

5-year trend of 

sulfate concentration, 

in microequivalents 

per liter

10-year trend of 

sulfate concentration, 

in microequivalents 

per liter

14-year trend of 

sulfate concentration, 

in microequivalents 

per liter

20-year trend of 

sulfate concentration, 

in microequivalents 

per liter

–4

–2

0

2

4

–2

2

–1

1

0

–0.5

0.5

1.5

–0.5

0

1.0

2.0

0

1.0

2.0

0.5

1.5

5.5

6.0

6.5

7.0

–0.05

–0.03

–0.01

–0.02

–0.04

0

0.01

–0.015

0

0.015

0.030

–0.005

0.010

0.025

0.040

–0.004

0

0.004

0.008

0.012

30

60

90

120

–1.5

0

1.5

3.0

–0.6

–0.2

–0.4

0.4

0

0.2

–1.5

–1.0

–0.5

0

0.5

–1.0

–0.6

–0.2

0.2

0.6

–6 –4 –2 2 4 6

Principal component 1

0 –4 –2 0 2 4 4

Principal component 2

–4 –3 –2 –1 0 1 2 3

Principal component 3

Figure 16. Scatterplot of principal component axes against water-quality variables, with 

significant correlations highlighted.

Figure 16. Principal component axes against water-quality variables, with significant 

correlations highlighted.
Analysis of Aquatic Macroinvertebrate Data 27

Analysis of Aquatic Macroinvertebrate 

Data

Park staff began sampling AM in 1984 with the objective 

of using measures of community structure and composition as 

indicators of stream condition (J.E.B. Wofford, National Park 

Service, written commun., 2011). AM monitoring is utilized 

by the park because these organisms, which are numerous and 

directly affected by changes in water-quality, are relatively 

easy to sample and interpret (National Park Service, 2008). 

The AM dataset for the park was used to evaluate the status 

and trends of AM communities within the study area. Geology 

and watershed size were examined as potential environmental 

predictors of AM communities because of their demonstrated 

importance for acidic deposition, which was also demonstrated 

by Cosby and others (2006), and macroinvertebrate community composition (Snyder and others, 2013).

Aquatic Macroinvertebrate Data Structure 

and Analysis

Park personnel have collected AM data since 1984; however, the current analysis was limited to data collected after 

1990 in order to facilitate integration of macroinvertebrate 

patterns with the physical and chemical assessments included 

herein. Furthermore, the analyses were limited to quantitative samples collected during the spring season (April–June) 

because summer sampling was discontinued in 1997, whereas 

spring sampling continued. Aquatic macroinvertebrate data 

derived from sampling that occurred in the spring of 1995 

were used to assess spatial patterns because more sites were 

sampled that year (N=54) than any other, and because limiting 

the analysis to a single year would allow control of temporal 

trends in this spatial analysis. The additional sampling effort 

that occurred in 1995 was largely the result of an increased 

effort—in the Paine Run, Staunton River, and Piney River 

watersheds—made with the intent of incorporating intensive 

sampling within representative watersheds of the park. This 

approach biased the results somewhat toward those watersheds, but allowed a better representation of geology and the 

gradient in watershed area.

Temporal trends were assessed using data collected at 

24 sites where approximately annual sampling occurred over 

the 20-year period between 1990 and 2009 (focal sites of 

LTEM program). Between 10 and 18 years of data were collected at focal sites over the 20-year period (table 8). Eleven of 

the 24 focal sites also were represented in the spatial dataset. 

Sites used in both spatial and temporal analyses are described 

in table 9 and are shown in figure 17. Because these sites are 

well distributed across gradients in watershed size and geology, inferences regarding landscape effects should be robust 

(fig. 18).

Stream sites were sampled using a portable invertebrate 

box sampler (PIBS) equipped with a 350-micrometer collecting net. At each site, three samples were collected at random 

locations within riffle areas, preserved in the field, and shipped 

to the laboratory of Dr. J. Reese Voshell at Virginia Polytechnic Institute and State University for taxonomic identification. 

Qualitative samples also were routinely collected, but not 

included in these analyses in order to ensure comparable levels 

of effort among samples, sites, and years for this assessment. 

Sampling protocols have remained consistent over the 20-year 

period of data collection, including personnel involved in 

taxonomic identification, allowing for robust assessment of 

spatial and temporal trends (J.E.B. Wofford, National Park 

Service, written commun., 2011). 

In the laboratory, individual specimens were identified 

to the “lowest practical level,” which was the species level 

of taxonomic resolution for some insect taxa, but the genus 

level for most insects, and higher levels (family or class) for 

most non-insect taxa, including worms and mollusks (J.E.B. 

Wofford, National Park Service, written commun., 2011). The 

variability in taxonomic resolution resulted in a large degree 

of taxonomic ambiguity in the raw database. That is, the same 

taxon often was identified at multiple levels of taxonomic 

resolution throughout the database, because of variation in 

the condition and maturity of individual specimens. This 

taxonomic ambiguity is common to large AM databases and 

can introduce additional uncertainty in computed summary 

measures (metrics) (Cuffney and others, 2007). 

For the analyses presented here, taxonomic ambiguity 

was addressed using a sequential process. First, all specieslevel identifications were lumped up to the genus level, as 

there was extensive ambiguity at the species level. Second, 

the “distribute parents among children” (DPAC) method in the 

Invertebrate Data Analysis System software (IDAS; Cuffney 

and Brightbill, 2011) was used to resolve the remaining ambiguity within a sample. The DPAC method assigns “ambiguous parents” (for example, specimens of the family Baetidae) 

to “children” (for example, genera Baetis and Pseudocleon) 

based on the proportions of the “children” observed in the 

sample. Resolving ambiguity within a sample ensured that 

assumptions were not made about the identity of ambiguous 

taxa that could compromise assessment of spatial and temporal trends.

Spatial and temporal trends in 14 AM metrics (table 10) 

were assessed. These metrics commonly are used either individually or in combination to infer stream condition (Rosenberg and Resh, 1993; Barbour and others, 1999). Metrics are 

largely based on richness, relative abundance, and ecological characteristics of individual taxa or groups of taxa. Taxa 

collected from park samples, and their associated ecological 

characteristics required for metric computation, are shown in 

appendix 1. PCA was used to assess correlations among AM 

metrics and was conducted separately for spatial and temporal 

datasets to determine correlations among AM metrics between 
28 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

EXPLANATION

!

!

Aquatic macroinvertebrate 

 study sites

Temporal

Spatial

Bedrock geology

Siliciclastic

Granitic

Basaltic

Shenandoah National 

 Park boundary

0 5 10 MILES

0 5 10 KILOMETERS

!

! !

!

!!

!(

!(

!(

!(

!

! ! ! !

!

!(

!(

!(

!( !(

!(

!(

!(

!(

!(

!(

!(

!(

!!

!(

!(

!(

!(

!(

!(

!

! !

!

!

!

!!

!!!!!!!!

!!

!

!

!

!!

!

!(

N

North District

South District

Central District

Figure 17. Map showing macroinvertebrate sampling sites in Shenandoah National Park. Location of study 

area provided in figure 1.

Figure 17. Macroinvertebrate sampling sites in Shenandoah National Park. Location of study area provided in 

figure 1.
Analysis of Aquatic Macroinvertebrate Data 29

sites within a year (spatial), as well as correlations among AM 

metric trend slopes (temporal). Prior to conducting PCA, the 

AM metrics were standardized using Z-scores so that differences in measurement scale among metrics would not skew 

PCA results (McCune and Grace, 2002). 

For analysis of spatial patterns, general linear modeling 

(GLM) was used to assess the influence of landscape setting, determined by contributing watershed area and bedrock geology, on macroinvertebrate metric means. Geology 

(siliciclastic, granitic, and basaltic) was a class variable in 

the analyses, indicated by the dominant geology upstream of 

the sample site, and watershed area was a continuous variable; the interaction between geology and watershed area also 

was evaluated. An interactive backward selection process 

for model building was used, beginning with the full model, 

which considered geologic effects, watershed area and their 

interaction. Non-significant terms were then removed one at 

a time, beginning with the interaction term. For interaction 

terms that were found to be significant predictors of a metric 

(p<0.10), the interaction was graphically evaluated to ensure 

the relation was not unduly influenced by outliers or points 

with high leverage. After the interaction term was evaluated, 

the same procedure was used to evaluate main effects. Finally, 

the relative importance of the predictor variables in the final 

model was evaluated using partial R2 statistics, which indicate 

the amount of variation in a metric explained by geology and 

watershed area. 

Temporal trends were assessed for four distinct 

time periods consistent with those used for water-quality 

trend analysis. Analysis of temporal trends included 

three steps. First, linear regression was used to determine 

site-specific relations (regression slopes, β) between each 

of the 14 AM metrics and year for each of the four time 

periods. For metrics that were not normally distributed, 

the data were transformed prior to performing linear 

regressions. All AM metrics that were not normally 

distributed were proportion or percent metrics; therefore, 

the arcsine square-root transformation was used, because 

it often improves normality of proportion data (Sokal and 

Rohlf, 1995). No attempt was made to evaluate nonlinear or threshold models in these exploratory analyses. 

Second, park-wide trends in each metric for each 

time period were evaluated using one-way t-tests. This 

analysis tested whether the mean trend slope was significantly different from zero (Sokal and Rohlf, 1995). 

Third, the influence of geology and watershed 

area on temporal trends was assessed. This was accomplished using essentially the same interactive backward 

model selection procedure used to assess spatial patterns (described earlier), except that trend slopes were 

the response variables in the analysis. For metric-time 

period cases in which the best model included geology as 

a significant predictor of temporal trends, a one-sample 

t-test was conducted within each geology class to test 

whether the mean trend slope was significantly different 

than zero. For metric time-period cases where the best 

model included watershed area as a significant predictor, linear 

regressions of trend slopes on watershed area were conducted 

to determine the percentage of variation in trend slopes (R2) 

explained by watershed area. Finally, for metric time-period 

cases in which the best model included the interaction between 

watershed area and geology as a significant predictor, linear 

regressions of trend slopes on watershed area were performed 

within each geology class.

Aquatic Macroinvertebrate Community 

Characteristics

The park’s AM dataset used to assess spatial patterns 

among 54 sites in 1995 contained records of 53,665 individual 

AMs and 105 unique taxa. The park’s AM dataset used to 

assess temporal trends among 24 sites between 1990 and 2009 

contained records of 500,781 individual macroinvertebrates 

represented by 140 unique taxa. In terms of taxonomic composition, the abundance of individuals in the orders Ephemeroptera, Plecoptera, and Trichoptera, broadly considered the three 

most generally sensitive taxonomic orders, composed about 

half of the total abundance in the database, although relative abundance varied between about 40 and 65 percent over 

the years (fig. 19A). Relative abundance of individuals in the 

family Chironomidae, generally considered a group relatively 

tolerant to stress, represented, on average, about one-third of 

Temporal Spatial Temporal Spatial Temporal Spatial

Log watershed area, in hectares 

1.5

2.0

2.5

3.0

3.5

Siliciclastic Granitic Basaltic

Geologic class

Figure 18. Graph showing watershed area of macroinvertebrate 

sampling sites used for temporal and spatial analyses, by geologic 

class.

Figure 18. Watershed area of macroinvertebrate sampling sites 

used for temporal and spatial analyses, by geologic class. 
30 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Table 8. Years in which samples were collected for 24 Long Term Ecological Monitoring sites used in temporal trend 

analysis. Site descriptions are provided in table 9.

[X, sample collected during year; —, sample not collected during year]

Site 

identifier

Years

of

data

Year

19901991199219931994199519961997199819992000200120022003200420052006200720082009

3F025 11 — — — — — — — X X X X X X X X — X — X X

3F103 16 X X X X X — — X X X X X X X X — X — X X

3F105 10 — — — — — — — — X X X X X X X — X — X X

3F109 10 — — — — — — — — X X X X X X X — X — X X

3F123 18 X X X X X X X X X X X X X X X — X — X X

3F143 11 — — — — — — — — X X X X X X X — X — X X

3F271 10 — — — — — — — — X X X X X X X — X — X X

3F300 18 X X X X X X X X X X X X X X X — X — X X

3F302 16 X X X X X — — X X X X X X X X — X — X X

1F002 17 X X X X X X — X X X X X X X X — X — X X

1F003 16 X X X X — — X X X X X X X X X — X — X X

2F009 16 X X X X X X — X — X X X X X X — X — X X

2F015 16 X X X X X X — X — X X X X X X — X — X X

2F072 18 X X X X X X X X X X X X X X X — X — X X

2F090 15 X X X X X — — X — X X X X X X — X — X X

2F131 18 X X X X X X X X X X X X X X X — X — X X

2F303 12 — — — — X X — X — X X X X X X — X — X X

2F306 18 X X X X X X X X X X X X X X X — X — X X

1F118 15 X X X X X — — X — X X X X X X — X — X X

1F308 18 X X X X X X X X X X X X X X X — X — X X

2F055 15 X X X X X — — X — X X X X X X — X — X X

2F253 12 — — — — — X X — X X X X X X X — X — X X

3F044 11 — — — — — — X — X X X X X X X — X — X X

3F084 11 — — — — — — X — X X X X X X X — X — X X
Analysis of Aquatic Macroinvertebrate Data 31

Table 9. Site information, watershed area, predominant geology, and site elevation for sites used in analysis of 

temporal trends and spatial patterns in aquatic macroinvertebrates.—Continued 

[ha, hectare; m, meter]

Site 

identifier Site name Watershedarea (ha)

Predominant

geology

Elevation

(m)

Analysis 

performed

1F003 Piney River 1,225 Granitic 366 Temporal

1F118 Jeremy’s Run 564 Basaltic 554 Temporal

2F055 Rose River & Hogcamp Branch 193 Basaltic 800 Temporal

2F090 Hazel River & tributary 275 Granitic 667 Temporal

3F025 Madison Run 1,420 Siliciclastic 430 Temporal

3F044 Moormans River, North Fork 1,240 Basaltic 498 Temporal

3F084 Moormans River, North Fork 2,522 Basaltic 341 Temporal

3F103 Twomile Run 561 Siliciclastic 370 Temporal

3F105 Lower Lewis Run 262 Siliciclastic 444 Temporal

3F109 Meadow Run 657 Siliciclastic 511 Temporal

3F271 White Oak Run 512 Siliciclastic 447 Temporal

3F302 Twomile Run 143 Siliciclastic 506 Temporal

2F072 Staunton River 1,065 Granitic 316 Temporal and spatial

1F002 Lands Run 171 Granitic 517 Temporal and spatial

1F308 Piney River - Upper 72 Basaltic 764 Temporal and spatial

2F009 Whiteoak Canyon Run 1,405 Granitic 343 Temporal and spatial

2F015 Rose River & Hogcamp Branch 2,368 Granitic 341 Temporal and spatial

2F131 Dry Run, N.F. 220 Granitic 493 Temporal and spatial

2F253 White Oak Canyon Run 234 Basaltic 950 Temporal and spatial

2F303 Hazel River 1,333 Granitic 329 Temporal and spatial

2F306 Staunton River 128 Granitic 900 Temporal and spatial

3F123 Paine Run 1,270 Siliciclastic 425 Temporal and spatial

3F143 Sawmill Run 618 Siliciclastic 463 Temporal and spatial

3F300 Paine Run 177 Siliciclastic 557 Temporal and spatial

1F001 Land’s Run 351 Granitic 360 Spatial

1F005 Piney River 986 Granitic 409 Spatial

1F006 Jeremy’s Run 2,202 Basaltic 285 Spatial

1F030 Thornton River, North Fork 1,832 Granitic 354 Spatial

1F137 Piney River 805 Basaltic 453 Spatial

1F138 Piney River 575 Basaltic 576 Spatial

1F316 Thornton River, South Fork 348 Granitic 436 Spatial

1F317 Overall Run 1,043 Siliciclastic 213 Spatial

2F035 Thornton River, South Fork 1,658 Granitic 292 Spatial

2F073 Staunton River 776 Granitic 412 Spatial

2F074 Staunton River 732 Granitic 462 Spatial

2F075 Staunton River 511 Granitic 619 Spatial

2F076 Staunton River 245 Granitic 753 Spatial

2F096 Pass Run 768 Basaltic 408 Spatial

2F252 Brokenback Run 1,005 Granitic 338 Spatial

2F315 Pass Run 1,076 Basaltic 335 Spatial

3F124 Paine Run 933 Siliciclastic 447 Spatial

3F125 Paine Run 793 Siliciclastic 458 Spatial

NPI02 Piney River 390 Basaltic 673 Spatial

NPR1 Paine Run 830 Siliciclastic 458 Spatial

NPR3 Paine Run 467 Siliciclastic 483 Spatial
32 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Table 9. Site information, watershed area, predominant geology, and site elevation for sites used in analysis of 

temporal trends and spatial patterns in aquatic macroinvertebrates.—Continued 

[ha, hectare; m, meter]

Site 

identifier Site name Watershedarea (ha)

Predominant

geology

Elevation

(m)

Analysis 

performed

PI02 Piney River 93 Granitic 393 Spatial

PI26 Piney River 145 Basaltic 779 Spatial

PI29 Piney River 65 Basaltic 692 Spatial

PR04 Paine Run 1,046 Siliciclastic 438 Spatial

PR05 Paine Run 93 Siliciclastic 450 Spatial

PR09 Paine Run 77 Siliciclastic 463 Spatial

PR15 Paine Run 737 Siliciclastic 471 Spatial

PR18 Paine Run 93 Siliciclastic 485 Spatial

PR19 Paine Run 596 Siliciclastic 483 Spatial

PR21 Paine Run 125 Siliciclastic 487 Spatial

PR25 Paine Run 466 Siliciclastic 485 Spatial

PR26 Paine Run 386 Siliciclastic 507 Spatial

PR28 Paine Run 351 Siliciclastic 509 Spatial

PR29 Paine Run 262 Siliciclastic 534 Spatial

PR31 Paine Run 260 Siliciclastic 535 Spatial

PR37 Paine Run 200 Siliciclastic 439 Spatial

SR04 Staunton River 181 Granitic 355 Spatial

SR09 Staunton River 624 Granitic 557 Spatial

SR11 Staunton River 102 Granitic 652 Spatial

SR13 Staunton River 49 Granitic 670 Spatial

Table 10. Description of the 14 macroinvertebrate community metrics evaluated.

[% , percent; <, less than; PTV, pollution tolerance value; SCI, stream condition index] 

Metric Definition Stress 

response

Richness Total number of taxa Decrease

EPT-r Total number of taxa in orders Ephemeroptera, Plecoptera, and Trichoptera Decrease

EPT%1 Percent abundance of individuals of orders Ephemeroptera, Plecoptera, and Trichoptera Decrease

E%1 Percent abundance of individuals in order Ephemeroptera Decrease

Hydro:T%1 Percent abundance of caddis flies (Trichoptera) in family Hydropsychidae (Trichoptera) Increase

Leuctra:P%1 Percent abundance of stoneflies (Plecoptera) represented by low pH tolerant genus Leuctra. Note, expected 

increase with acid stress only

Increase

DOM2 Relative abundance (%) of the two most abundant taxa Increase

SimpsonD Simpson diversity index (1-D). Measure of richness weighted by evenness Decrease

PTV Pollution tolerance value. Weighted sum of total taxa by taxon-specific tolerance values (Klemm and others, 2002) Increase

Intol%1 Percent abundance of macroinvertebrates with tolerance values <2 Decrease

Scraper%1 Percent abundance of functional feeding group containing “scrapers” (Merritt and Cummins, 1996) Decrease

Shredder%1 Percent abundance of functional feeding group containing “shredders” (Merritt and Cummins, 1996) Decrease

Hapto%1 Haptobenthos. Percent abundance of macroinvertebrates requiring clean, firm, coarse substrates (“crawlers” + 

“clingers”) (Merritt and Cummins, 1996)

Decrease

SCI1 Stream condition index. Multi-metric index comprising 8 individual metrics developed for highland streams in 

Virginia (Burton and Gerritsen, 2003)

Decrease

1

Metric was arcsine square-root transformed prior to analyses.
Analysis of Aquatic Macroinvertebrate Data 33

Ephemeroptera

Plecoptera

Trichoptera

Chironomidae

Other insects

Non-insects

Taxonomic category

EXPLANATION

Predator

Collector-filterer

Collector-gatherer

Shredder

Scraper

Trophic category

EXPLANATION

199019911992199319941995199619971998

1999

2000200120022003

20042005

2006

2007

20082009

Year

0

20

40

60

80

100

B. Trophic composition

Percent abundance

0

20

40

60

80

100

Figure 19. Graphs showing A, taxonomic, and B, trophic composition of 

aquatic macroinvertebrates by sample year computed from samples 

collected from 24 sites.

A. Taxonomic composition

Figure 19. A, Taxonomic and B, trophic composition of aquatic macroinvertebrates by sample 

year computed from samples collected from 24 sites.
34 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

the total AM abundance. Although it is difficult to make direct 

comparisons between patterns of taxonomic composition 

observed in the park and other datasets from around the region 

because of differences in sampling methods, the relative abundances of sensitive and tolerant taxonomic groups observed in 

the park suggest that, overall, streams in the park are representative of high-quality headwater streams in the region. Other 

insect taxa (orders Hemiptera, Megaloptera, Odonata, and 

non-midge Diptera) also represented a significant taxonomic 

component (mean=15 percent), whereas non-insect taxa were 

a relatively minor component of the AM fauna in terms of 

abundance (fig. 19A). 

Stream communities are structured by the source and 

quality of food and nutrient inputs entering streams. These 

inputs are linked to characteristics of the stream and contributing watershed, such as stream width and riparian shading, 

that change in predictable ways as streams progress from 

headwaters to mouth (Vannote and others, 1980). Metrics that 

describe trophic composition are often included in stream condition assessments because they confer information regarding 

food web structure. That structure may be disrupted by natural 

and anthropogenic disturbances, including watershed deforestation, pest-induced changes in forest composition, and nutrient enrichment through land use or atmospheric deposition 

(Likens, 2010). In terms of trophic composition, collectors 

(comprising collector-filterer and collector-gatherer trophic 

groups) were the dominant trophic group, representing slightly 

more than half of the total number of individuals collected 

over the 20 years at the 24 sites (fig. 19B). Of the total annual 

AM abundance, (on average) shredders composed about 

17 percent, scrapers 11 percent, and predators about 9 percent 

(fig. 19B). This distribution of trophic groups is fairly typical 

for larger (2nd and 3rd order) headwater streams (Haggerty 

and others, 2002; Heino 2005).

PCA results indicated that there were strong correlations 

among AM metric site scores among sites (spatial dataset; 

table 11), and among AM metric trend slopes (temporal dataset; table 11). Over 83 percent of the within-year spatial variation was explained by three PCA gradients, and nearly 77 percent of the temporal variation in AM metrics was explained 

by three PCA gradients (table 11). High correlations among 

the 14 metrics evaluated are not surprising, given that some 

groups of metrics are proportional and related (for example, 

Shredder% and Scraper%), and other metrics are simply subsets (for example, Richness and EPT-r). Distinct differences in 

the identity of correlated metrics between the spatial and temporal datasets, however, indicate that the factors that structure 

spatial variation are different than the factors that drive temporal trends. For example, spatial variation in the functional 

metrics Shredder% and Scraper% were highly correlated with 

the integrated stream condition index (SCI) and E% (table 11), 

whereas temporal changes in those metrics were uncorrelated 

(table 11). For this reason, spatial and temporal trend analyses 

were conducted on the full complement of metrics. 

Table 11. Factor loadings from principal components analysis of aquatic macroinvertebrate metrics for data used in the spatial and 

temporal analyses of aquatic macroinvertebrate communities.

[ns, correlation not significant for α=0.10; SCI, stream condition index; PTV, pollution tolerance value]

Spatial dataset Temporal dataset

Metric Principal 

component 1

Principal 

component 2

Principal 

component 3 Metric Principal component 1

Principal 

component 2

Principal 

component 3

SCI 0.936 ns ns PTV –0.943 ns ns

DOM2 –0.924 ns ns Intol% 0.92 ns ns

SimpsonD 0.917 ns ns EPT% 0.911 ns ns

Leuctra:P% –0.847 ns ns Hapto% 0.842 ns ns

E% 0.834 ns ns SCI 0.825 ns ns

Scraper% 0.812 ns ns SimpsonD 0.791 ns ns

Shredder% –0.748 ns ns DOM2 –0.707 ns ns

Hapto% ns 0.778 ns E% 0.642 ns ns

PTV ns –0.938 ns Shredder% ns –0.93 ns

Intol% ns 0.906 ns Leuctra:P% ns –0.752 ns

EPT% ns 0.865 ns Scraper% ns 0.532 ns

Richness ns ns 0.93 Richness ns ns 0.928

EPT-r ns ns 0.915 EPT-r ns ns 0.884

Hydro:T% ns ns ns Hydro:T% ns ns –0.533

Variance explained 40.51 26.36 16.17 Variance explained 41.84 19.07 16.02
Analysis of Aquatic Macroinvertebrate Data 35

Spatial Patterns in Aquatic Macroinvertebrates

A large fraction of spatial variation in AM metrics 

derived from samples collected in 1995 was explained by 

geology and watershed area (table 12). Of the 14 metrics 

assessed, 13 showed significant relations with one or both 

landscape predictors. Only the metric Hydro:T% showed no 

relation with either watershed area or geology. The amount 

of variation explained by these two variables was substantial, 

ranging from weak relations such as EPT%, for which less 

than 15 percent of total variation was explained, to very strong 

relations such as Shredder%, for which nearly 80 percent of 

the total spatial variation was explained by geology and watershed area (table 12).

In general, geology was a stronger predictor than watershed area for most standard metrics (those not excluding 

Leuctra). Geology was a significant predictor of 13 metrics, 

with 6 metrics showing significant relations with geology 

alone, 6 metrics showing significant but independent (that is, 

no interaction between geology and watershed area) relations 

with geology and watershed area, and 1 metric significantly 

influenced by the interaction between geology and watershed area (table 12). The amount of variation in AM metrics 

explained by geology ranged from 14 percent for Intol% to 

64 percent for SimpsonD (table 12). In contrast, watershed 

area was a significant predictor for only seven AM metrics, 

and in only one case (Intol%) was it found to be more important than geology (table 12). 

The observation that geology was a stronger predictor of 

AM community metrics is in contrast to what Snyder and others (2013) found when evaluating the importance of headwater 

streams in the park, where watershed area was found to be 

more important than geology for predicting the same set of 

metrics. In that assessment, however, small headwater streams 

(<100-hectare drainage areas) and perennial springs were 

included, whereas this analysis included long-term AM sampling in a small fraction of sites with watershed areas less than 

100 hectares (ha). These differences between the two studies 

suggest a non-linear, bimodal relation between AM metrics and 

watershed area, wherein relatively large changes in AM communities occur as the contributing watershed area decreases 

to less than 100 ha. Snyder and others (2013) hypothesized 

that the reduced stream condition inferred from AM metrics 

in headwater areas was due to the increased likelihood of flow 

disturbance (especially droughts) in smaller watersheds.

The observed relations between geology and AM assemblage structure suggests a strong influence of water quality on 

AM assemblage structure, because ANC and pH have shown to 

be strongly tied to watershed geology in the park by this report 

and others (Herlihy and others, 1993; Cosby and others, 2006; 

Snyder and others, 2006). In general, the relations between 

most AM metrics and geology suggest that stream condition 

was poorer in siliciclastic watersheds where acidity was highest because of the low buffering capacity of rocks and soils 

in this geology type (Cosby and others, 2006). For example, 

most macroinvertebrate metrics expected to decline with stress 

(Richness, EPT-r, SimpsonD, E%, Hapto, SCI, and Scraper%) 

were lower in siliciclastic watersheds, and metrics expected to 

increase with stress (DOM2) were higher in siliciclastic watersheds (fig. 20). In addition, the Leuctra:P% metric, which was 

adopted by the park’s managers because of the known tolerance of the stonefly genus Leuctra to acidity (Kimmel and others, 1985; Fjellheim and Raddum, 1990), also showed a pattern 

that suggests acid stress was most pronounced in siliciclastic 

watersheds. Specifically, a greater proportion of Leuctra was 

present in siliciclastic watersheds than in others (fig. 20). 

In contrast to the prevailing pattern, three AM metrics—

EPT%, Intol%, and pollution tolerance value (PTV)—showed 

the opposite pattern in metric scores, thereby suggesting that 

stream condition was better at sites in siliciclastic watersheds 

than in others. Both EPT% and Intol%, metrics generally 

expected to decline with stress, were highest in siliciclastic 

watersheds; moreover, PTV, a metric typically expected to 

increase with stress, was lowest in siliciclastic watersheds 

(fig. 20). These three metrics deviate from the prevailing 

pattern because of the overwhelming influence of a single 

taxon, the stonefly (Leuctra sp) and its relative sensitivity to 

different stressors. 

Table 12. Partial coefficients of determination (partial R2) for 

variables found to be statistically significant predictors (p<0.05) 

for each metric, and overall model R2, from interactive stepwise 

general linear modeling used to assess the effects of geology and 

watershed area on spatial variation in macroinvertebrate metrics.

[NI, not interpretable (main effects not interpretable when interactions are significant); ns, not significant; R2

, coefficient of determination; PTV, pollution 

tolerance value; SCI, stream condition index]

 Metric

Partial R2

Model 

R2 Area Geology Area × 

geology

Richness 0.159 0.316 ns 0.475

EPT-r 0.179 0.316 ns 0.495

EPT% ns 0.145 ns 0.145

EPT% (without Leuctra) ns 0.402 ns 0.402

E% ns 0.445 ns 0.445

Hydro:T% ns ns ns 0

Leuctra:P% ns 0.318 ns 0.318

DOM2 ns 0.636 ns 0.636

SimpsonD ns 0.64 ns 0.64

PTV 0.115 0.207 ns 0.322

PTV (without Leuctra) 0.153 0.229 ns 0.382

Intol% 0.196 0.136 ns 0.332

Intol% (without Leuctra) 0.146 0.382 ns 0.528

Scraper% ns 0.579 ns 0.579

Shredder% NI NI 0.138 0.783

Shredder% (without Leuctra) 0.200 ns ns 0.200

Hapto% 0.121 0.172 ns 0.293

SCI 0.091 0.583 ns 0.674
With Leuctra EXPLANATION Without Leuctra With Leuctra EXPLANATION Without Leuctra With Leuctra EXPLANATION Without Leuctra SimpsonD 0.6 0.7 0.8 0.9 1.0 Leuctra:P% 0 0.2 0.4 0.6 0.8 1.0 DOM2 20406080 100 Richness 102030405060 EPT-r 10152025303540

Siliciclastic Granitic Basaltic

PTV

123456

Siliciclastic Granitic Basaltic

Intol%

0

0.20.40.60.8

Hapto%

0.40.50.60.70.80.9

Siliciclastic 

SCI

0.50.60.70.80.9

Scraper%

0

0.10.20.30.40.50.6

E%

0

15304560

A B B A B B A C B C 

A B B B A B

A B B B A B B A A A B A

A A 

A C B

Siliciclastic Granitic

Geologic class Geologic class Geologic class 

Basaltic

EPT%

0.10.20.30.40.50.60.70.8

B A B B A B

Figure 20. Graphs showing aquatic macroinvertebrate metric values by geologic class for 12 metrics significantly influenced by geology. Legraph refer to the results of Tukey’s Multiple range tests for differences in group means (means for groups with the same letter are not signieach other at α=0.05). The raw (untransformed) data were plotted but statistical modeling was conducted on transformed data for metrics w

Aquatic macroinvertebrate metric

Figure 20. Aquatic macroinvertebrate metric values by geologic class for 12 metrics significantly influenced by geology. Letters at the refer to the results of Tukey’s multiple range tests for differences in group means (means for groups with the same letter are not significeach other at α=0.05). The raw (untransformed) data were plotted but statistical modeling was conducted on transformed data for metridistributions.
Analysis of Aquatic Macroinvertebrate Data 37

Although Leuctra is tolerant to acidic conditions (Kimmel and others, 1985; Fjellheim and Raddum, 1990), it is 

widely considered to be highly sensitive to anthropogenic 

stress in general—the tolerance value for Leuctra is zero 

(appendix 1), which is the most sensitive assignment possible 

on the 0–10 scale. Additionally, Leuctra is one of the most 

stenothermic genera of stonefly taxa, requiring a relatively 

narrow range of cool or cold water to complete their life cycle, 

and has been implicated as particularly vulnerable to warming water temperatures (Sweeney and others, 1992; Fochetti 

and Tierno De Figueroa, 2008; Tierno de Figueroa and others, 

2010). Tolerance values were originally derived and modified 

from the seminal work by Lenat (1993) that based tolerance 

values on taxon-specific abundances across a water-quality 

gradient in North Carolina streams. These values, with 

regional modifications, are widely used in AM multimetric indices used to infer stream condition throughout North 

America (Barbour and others, 1999; Klemm and others, 2002; 

Burton and Gerritsen, 2003) and have been adopted by the 

park’s management. Although Lenat (1993) does not describe 

the types of water-quality issues represented in the North 

Carolina state database, it is likely that the predominant stressors were sediments and chemical pollutants associated with 

industrial and municipal wastewater as well as runoff from 

urban and agricultural development. High-elevation watersheds vulnerable to acidic deposition were probably not well 

represented, and consequently, it is unlikely that acidity was 

a major water-quality issue in the streams used to compute 

tolerance values. In this study, taxon-specific tolerance values 

were used to compute the metric PTV directly and identify 

intolerant taxa for computing the Intol% metric (see table 10). 

Therefore, in regions where stream acidification is an important issue, it may be inappropriate and counterproductive to 

include the acid-tolerant taxa Leuctra in these AM metrics, 

which are designed to be indicative of general stream condition, especially if Leuctra is a numerically abundant taxon. 

In the park, Leuctra was the second most abundant taxon 

collected. More than 55,000 Leuctra individuals out of a 

nearly 501,000 total specimens (>10 percent) were collected at 

the 24 LTEM sites over the 20-year monitoring period. Only 

individuals of the family Chironomidae were more abundant, 

and the family Chironomidae comprises more than a dozen 

individual genera in the park. Moreover, Leuctra composed a 

significant fraction of the EPT and intolerant taxa assemblages 

in the park, and the acid-tolerant genus was more abundant in 

siliciclastic watersheds than in others (fig. 21). Thus, given its 

EXPLANATION

90th percentile

75th percentile

Median

25th percentile

10th percentile

Outlier

Outlier

Community

Siliciclastic Granitic Basaltic

% Leuctra 

0

20

40

60

80

100

EPT taxa

Siliciclastic

Geologic class Geologic class Geologic class Geologic class

Granitic Basaltic

Intolerant taxa

Siliciclastic Granitic Basaltic

Shredder taxa

Siliciclastic Granitic Basaltic

Figure 21. Percent abundance of Leuctra relative to the overall community, the Ephemeroptera, Plecoptera, and 

Trichoptera (EPT) assemblage, the intolerant taxa assemblage, and the shredder assemblage (right) by watershed 

type. Boxplots report the median (black line within boxes), the 25th and 75th percentiles (limits of box), the 10th and 

90th percentiles (limits of error bars), and all outliers (points outside the 10th and 90th percentiles) as points.
38 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

numerical abundance, it is clear that considering the acidtolerant genus Leuctra, a sensitive taxa based on its tolerance 

value, has the potential to greatly confound interpretations of 

trends in ecological condition.

Based on these findings, Leuctra was removed from the 

calculations of EPT%, Intol%, and PTV metrics and spatial 

patterns were reevaluated. In all three cases, the strength of 

the relations between geology and AM metrics increased 

(table 12), indicating that the Leuctra-adjusted metrics more 

closely reflected geologically-induced variation in water 

quality. Moreover, after removing Leuctra, these metrics were 

consistent with the prevailing pattern observed in the other 

metrics. That is, metric scores indicated that stream condition was poorest in siliciclastic watersheds (fig. 20). Taken 

together, the spatial patterns in AM metric scores correspond 

well to geologically-induced gradients in acidity, and to previously reported relations between acidity and aquatic macroinvertebrates (Kimmel and others, 1985; Griffith and others, 

1995; Bowman and others, 2006). 

Observed spatial relations between AM assemblage structure and geology were not entirely consistent with an acidity 

gradient. Specifically, although AM metric scores in siliciclastic watersheds, which have the lowest mean pH and ANC, 

were always poorest, AM metric scores in basaltic watersheds 

tended to be worse than scores in granitic watersheds (fig. 20), 

even though mean pH and ANC are typically higher in basaltic 

than in granitic watersheds, as demonstrated in this study and 

by Cosby and others (2006). Spatial patterns in AM metrics 

correlate better with total SO4

2– concentrations than with pH 

and ANC. Mean SO4

2– concentrations were highest in siliciclastic watersheds, followed by basaltic and then granitic 

watersheds, an effect also documented by Cosby and others 

(2006) and Snyder and others (2013). Although high levels 

of sulfate salts are known to cause osmoregulation problems 

and toxicity in freshwater macroinvertebrates (Soucek, 2007), 

streams in the park are low-ionic strength waters having 

specific conductance values typically ranging between 18 and 

50 microsiemens per centimeter (μS/cm); consequently, SO4

2–

concentrations are not toxic at observed levels. 

Sulfate is the primary anion responsible for acidification 

in mid-Atlantic region streams (Herlihy and others, 1993). The 

observation that pH and ANC are higher in basaltic watersheds 

than granitic watersheds, despite greater SO4

2– concentrations, 

suggests that the concentrations of base cations (calcium, 

magnesium, sodium, and potassium) in basaltic watersheds 

are sufficient to neutralize the effect of the increased SO4

2–, 

at least on average. Based on these observations, it is hypothesized that, although average pH and ANC values are higher in 

basaltic watersheds, the higher SO4

2– concentrations in basaltic 

watersheds lead to acidification levels greater than in granitic 

watersheds during some biologically important times such as 

storm events. In this case, although mean ANC and pH might 

be lower in granitic watersheds, minimum pH and ANC would 

be expected to be lower in basaltic watersheds, and these 

most extreme conditions may be more important than average conditions to AM populations. A more detailed analysis of 

storm-event water-quality patterns would be required to test 

this hypothesis. Although data for such an analysis are available, the analysis was beyond the scope of the current effort.

In addition to geology, watershed area was found to be 

a significant predictor of spatial variation in six AM metrics, 

although the relation was relatively weak in all cases (partial 

R2<0.20, table 12). The two richness metrics (Richness and 

EPT-r) both increased with watershed area (fig. 22). This 

pattern is not surprising and generally conforms to patterns 

observed by others, whereby macroinvertebrate richness generally increases from headwaters to mid-reach sized streams 

(Clarke and others, 2008, and references therein). 

Spatial patterns in PTV, Intol%, Hapto%, and SCI suggest 

stream condition was better in smaller watersheds than larger 

watersheds. Specifically, PTV, a metric generally expected to 

increase with stream-condition degradation, increased with 

watershed area; in contrast, Intol%, Hapto%, and SCI, metrics 

expected to decline with stream-condition degradation, all 

decreased with watershed area (fig. 22). Unlike the case for 

watershed geology, the relation between the Intol% metric and 

watershed area remained essentially unchanged after removing 

Leuctra (table 12). Once again, this is the opposite of what Snyder and others (2013) found for these four metrics in their comparison of small headwater streams to larger streams. In that 

study, the authors found that patterns for these same four AM 

metrics suggested that stream condition improved as watershed 

area increased. This discrepancy suggests a non-linear and possibly bimodal relation between watershed area and many AM 

community metrics, whereby the relation changes as contributing watershed area declines to less than 100 ha. In addition, 

although all six of these metrics showed significant relations 

with geology (fig. 20), the fact that the interaction between 

watershed area and geology was not a significant predictor for 

any of them suggests that changes in AM metrics with watershed area were probably not associated with corresponding 

changes in water quality. Based on these factors, it is hypothesized that the weak relations between the four non-richness 

metrics of stream condition and watershed area observed from 

this analysis are probably coincidental, relating more to natural 

changes in component taxa habitat preferences with changing 

watershed area than to some form of degradation of habitat in 

larger watersheds. These results emphasize the need to account 

for natural variation in physical factors such as watershed area 

when interpreting bioassessment results. 

Spatial variation in one metric, Shredder%, was related to 

the interaction between watershed area and geology (table 12). 

Shredder% showed a highly significant decline with watershed area in siliciclastic watersheds, a weak but significant 

increase with watershed area in granitic watersheds, and no 

relation with watershed area in basaltic watersheds (fig. 23A). 

The relation observed in siliciclastic watersheds is consistent 

with what is typically expected in unperturbed watersheds, in 

that shredders are typically most abundant in headwater areas 

where streams are narrow, leaf inputs are high, and light levels 

are low, and gradually decline as streams widen (Vannote and 

others, 1980; Wohl and others, 1995). 
Analysis of Aquatic Macroinvertebrate Data 39

With Leuctra

EXPLANATION

Without Leuctra 

With Leuctra

EXPLANATION

Without Leuctra 

Richness

10

20

30

40

50

60

1.5 2.0 2.5 3.0 3.5

Hapto%

0.4

0.5

0.6

0.7

0.8

0.9

Pollution tolerance value (PTV)

1

2

3

4

5

6

EPT-r

5

10

15

20

25

30

35

1.5 2.0 2.5 3.0 3.5

Stream condition index (SCI)

0.55

0.60

0.65

0.70

0.75

0.80

0.85

Intol%

0

0.2

0.4

0.6

0.8

Log watershed area, in hectares Log watershed area, in hectares

Figure 22. Graphs showing relations between aquatic macroinvertebrate (AM) metrics and 

watershed area for the six metrics significantly influenced by watershed area.

Aquatic macroinvertebrate metric

Figure 22. Relations between aquatic macroinvertebrate (AM) metrics and watershed area for 

the six metrics significantly influenced by watershed area. 

The Shredder% metric was unduly influenced by the acidtolerant stonefly, Leuctra, as with PTV, Intol% and EPT% metrics discussed above. Leuctra was the dominant shredder taxon 

in the park, composing as much as 100 percent of the shredder 

assemblage in some samples. As with the other three metrics, 

the relative abundance of Leuctra in the shredder assemblage 

was greatest in siliciclastic watersheds (fig. 21). When Leuctra

was removed and the Shredder% metric recomputed, geology 

was no longer a significant predictor of Shredder% (table 12), 

and the relation between Shredder% and watershed area 

reflected the expected pattern of gradual decline in Shredder% 

with watershed area for all geology types (fig. 23B). 

Temporal Trends in Aquatic Macroinvertebrates

Temporal trends for all 14 AM metrics were observed, 

although trends were park-wide for some metrics, whereas 

they depended on underlying geology or watershed area for 

others (table 13). In addition, trends seemed to be relatively 

consistent among time periods for the seven metrics that 

exhibited park-wide trends (E%, Leuctra:P%, Hapto%, Richness, EPT-r, SimpsonD, and DOM2), although the magnitude 

and statistical significance of trends varied. In contrast, for 

those metrics whose trends depended on landscape factors, the 

relations often differed between time periods. For example, the 
40 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

1.5 2.0 2.5 3.0 3.5

0

0.2

0.4

0.6

1.5 2.0 2.5 3.0 3.5 1.5 2.0 2.5 3.0 3.5

Shredder%

metric score

A. With Leuctra 

1.5 2.0 2.5 3.0 3.5

0

0.2

0.4

0.6

1.5 2.0 2.5 3.0 3.5 1.5 2.0 2.5 3.0 3.5

B. Without Leuctra 

Siliciclastic class Granitic class

Siliciclastic class Granitic class

Basaltic class

Basaltic class

Watershed area, in hectares

Figure 23. Graphs of Shredder% computed A, with Leuctra, and B, without Leuctra, against watershed area, 

by geologic class.

Figure 23. Shredder% computed A, with Leuctra, and B, without Leuctra, against watershed area, by 

geologic class.

Table 13. Description of temporal trends of 14 aquatic macroinvertebrate 

metrics. 

[Statistically significant (p<0.10) temporal trends are identified as park-wide or dependent upon underlying geology or watershed area for each time period evaluated. ns, not 

significant; PTV, pollution tolerance value; SCI, stream condition index]

Metric 20-year 

trend

14-year 

trend

10-year 

trend

5-year 

trend

Richness ns Park-wide Park-wide Park-wide

Ept-r ns Park-wide Park-wide Park-wide

EPT%1 Geology ns ns Watershed area

E% Park-wide ns ns Park-wide

Hydro:T% Geology ns ns ns

Leuctra:P Park-wide Park-wide Park-wide Park-wide

DOM2 ns ns Park-wide ns

SimpsonD ns ns Park-wide ns

PTV1 Geology ns ns Watershed area

Intol%1 Geology Geology ns Watershed area

Scraper% Geology ns ns Park-wide

Shredder%1 ns Geology ns Watershed area

Hapto%1 Park-wide Park-wide Park-wide Park-wide

SCI Geology ns ns ns

1

The genus Leuctra was not included in calculating these metrics.
Analysis of Aquatic Macroinvertebrate Data 41

20-year trends for six AM metrics (EPT%, Hydro:T%, PTV, 

Intol%, Scraper%, SCI) were significantly influenced by geology (table 13), suggesting a relation to water-quality trends. In 

contrast, the more recent 5-year trends for three of those same 

six metrics (EPT%, PTV, Intol%) depended on watershed area 

and were not influenced by geology (table 13). These patterns suggest that the environmental factors driving long-term 

trends (20-year) differed from those responsible for short-term 

(5-year) trends for these metrics. 

Preliminary review of the AM metrics exhibiting parkwide trends indicated no clear temporal trends in terms of 

overall stream condition, with some metrics seeming to 

indicate improving stream condition and others suggesting 

declines. For example, AM diversity, as indicated by the Richness and EPT-r metrics, declined significantly over the 5-, 10-, 

and 14-year time periods, and the SimpD and DOM2 metrics 

(also measures of diversity or evenness) showed significant 

declines in stream condition, but only for the 10-year time 

period (fig. 24). Likewise, the Leuctra:P% metric showed 

highly significant declines for all four time periods evaluated 

(fig. 24). This metric could be indicative of either improving (declining) stream acidity, because Leuctra is relatively 

tolerant to acidic conditions, or declining stream conditions, 

because Leuctra is also known to be highly sensitive to other 

chemical pollutants and habitat degradation. 

In contrast to the metrics indicating declining conditions, 

Hapto% and E%, two AM metrics expected to decline with 

stress, showed significant increasing trends over time, suggesting improvement in stream condition. The Hapto% metric 

showed increasing trends for all four time periods evaluated, 

and the E% metric showed significant increases for both longterm (20-year) and short-term (5-year) time periods (fig. 24). 

After excluding the genus Leuctra, however, the Hapto% 

metric was mainly composed of mayflies (E%), thus the strong 

similarity in temporal patterns of the two metrics (fig. 24). In 

addition to the proportion of mayflies (E%), temporal trends in 

the actual abundance of mayflies were evaluated (not shown) 

and no significant trends were detected. Consequently, the 

increasing trends observed for both of these metrics were 

largely due to the decreasing abundance of other taxa (particularly the genus Leuctra); therefore, these metrics are probably 

not indicative of improving stream condition.

For other AM metrics, temporal trends were related to 

geology and watershed area. Geology was found to have 

a significant effect on temporal trends for seven AM metrics, although the influence was limited to long-term trends 

(fig. 25). Variation in long-term trends explained by geology 

ranged from about 20 percent for the Hydro:T% metric to 

over 50 percent for Intol% (fig. 25). In general, where geology 

was an important predictor of trends, metric scores at sites in 

siliciclastic watersheds indicated declining stream condition 

over time relative to sites in granitic and basaltic watersheds 

(EPT%, Hydro:T%, PTV, Intol%, Scraper%) or, as in the case 

of SCI, sites in granitic and basaltic watersheds indicated 

improving stream condition relative to sites in siliciclastic 

watersheds (fig. 25). The effect of geology was also evident 

at the 14-year timeframe for Intol%, although substantially 

less variation was explained than at the 20-year timeframe 

(fig. 25). Geologic effects on Shredder% at the 14-year timeframe were comparable to those for Intol% (fig. 25). Again, 

in both of these cases, sites in siliciclastic watersheds showed 

evidence of declining stream condition relative to sites in 

granitic and basaltic watersheds. 

Watershed area also affected trends for four AM metrics (fig. 26) although, in contrast to underlying geology, the 

influence of watershed area was limited to short-term (5-year) 

trends (not shown). Patterns for three of the four AM metrics 

(Shredder%, PTV, and Intol%) suggest that stream condition 

degraded between 2005 and 2009 in smaller watersheds (less 

than 1,000 ha) and either improved or did not change in larger 

watersheds (fig. 26). Similarly, the 5-year trend for the EPT% 

metric suggests that stream condition did not change significantly in smaller watersheds but improved in larger watersheds (fig. 26). 

In summary, the prevailing evidence derived from AM 

metric trends suggests an overall decline in stream condition over the last 20 years in the park, although a few metrics 

showed some evidence of improvement in stream condition in 

selected landscape settings. Specifically, the PTV, Intol%, and 

SCI metrics showed evidence of improving trends in streams 

underlain by granitic and (or) basaltic bedrock but declining 

trends in streams underlain by siliciclastic bedrock. Declining trends in macroinvertebrate community conditions were 

evident at both short-term and long-term time scales and, 

although some metrics showed declines park-wide, declines 

in other metrics depended on landscape factors (geology type, 

watershed area, or correlates). 

Additional information can be derived from the trend 

data to inform those designing new or additional AM trend 

monitoring efforts in the park’s streams about the period of 

data collection required to detect changes in the AM community. For a given metric, the trend slope alone (even if 

statistically significant) does not provide meaningful information regarding the relative sensitivity or statistical power of 

the metric to detect change over time. Conventional power 

analysis equations (Zar, 1999), therefore, were used to compute the minimum detectable change (MDC), which is defined 

as the absolute change in the mean metric value that can be 

detected between two samples (in this case, time periods), 

given observed variances. MDC cannot be directly compared 

among metrics, however, because metrics vary with respect 

to scale. The number of years required to detect a significant 

change (YTDC) was computed by dividing MDC by the mean 

trend slope. YTDC provides a sense of the magnitude of the 

observed trends and how they compare among metrics. 

For the aforementioned analyses, YTDC was computed 

for two sets of metrics. For the five metrics that exhibited 

park-wide trends, data from the 10-year timeframe were used 

because trend slopes, which indicate change in metric value 

per year, were greatest for the 2000–2009 time period. The 

Hapto% metric is not shown because the long-term trend for 

this metric was suspected to be spurious, as discussed earlier. 
42 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Hapto%

–0.05

0

0.05

0.10

0.15

0.20

Richness

–8

–6

–4

–2

0

2

EPT-r

–6

–4

–2

0

2

E%

–0.04

0

0.04

0.08

0.12

Leuctra:P%

–0.05

0

0.05

0.10

SimpsonD

–0.2

–0.1

0

0.1

0.2

Trend slope

Time period, in years

Time period, in years

DOM2

–12

–6

0

6

12

Figure 24. Graphs of trend slopes for the seven aquatic macroinvertebrate (AM) metrics with trends independent 

of geology and watershed area (park-wide) by time period. p-values from one-sample t-tests to determine 

significance (α = 0.10) are shown above the plotted points. Data used to compute trend slopes were arc sine 

transformed for the Leuctra:P%, Hapto%, and E% metrics; therefore, plotted slope values cannot be used to infer 

percent change per year directly. The trend slopes for the Richness, EPT-r, SimpsonD, and DOM2 metrics were 

computed from regressions of untransformed data and therefore can be interpreted directly.

20

0.004

20

0.812 0.216

20

0.003

20

0.000

20

0.326

20

0.678

20

14

0.046

14

0.068

14

0.007

14

0.111

14

0.000

14

0.323

14

0.513

10

0.001

10

0.000

10

0.000

10

0.528

10

0.000

10

0.059

10

0.040

5

5

5

5

0.001

5

0.000

5

0.171

5

0.003

0.000

0.000

0.147

Figure 24. Trend slopes for the seven aquatic macroinvertebrate (AM) metrics with trends 

independent of geology and watershed area (park-wide) by time period. p-values from 

one-sample t-tests to determine significance (α = 0.10) are shown above the plotted points. 

Data used to compute trend slopes were arc sine transformed for the Leuctra:P%, Hapto%, 

and E% metrics; therefore, plotted slope values cannot be used to infer percent change per 

year directly. The trend slopes for the Richness, EPT-r, SimpsonD, and DOM2 metrics were 

computed from regressions of untransformed data and therefore can be interpreted directly.
Analysis of Aquatic Macroinvertebrate Data 43

EPT%

–0.03

–0.02

–0.01

0

0.01

0.02

Hydro:T%

–0.02

–0.01

0

0.01

0.02

PTV

–0.1

0

0.1

0.2

Intol%

–0.02

–0.01

0

0.01

0.02

Scraper%

Siliciclastic Granitic Basaltic

–0.010

–0.005

0

0.005

0.010

0.015

SCI

Siliciclastic Granitic Basaltic

–0.010

–0.005

0

0.005

0.010

0.056 0.486 0.114 0.127 0.354 0.320

0.080 0.055 0.013 0.095 0.046 0.003

0.017 0.308 0.134

0.315 0.050 0.013

R2

 = 0.289

R2

 = 0.434

R2

 = 0.204

R2

 = 0.503

R2

 = 0.296

0.051 0.339 0.386

R2

 = 0.192

R2

 = 0.302

Shredder%

–0.02

–0.01

0

0.01

0.02

R2 = 0.304

0.045 0.550 0.208

Geologic class

Geologic class

Trend slopes

Figure 25. Graphs of trend slopes by geologic class for the seven metrics with long-term trend slopes indicating significant 

geologic effects. The p-values from one-sample t-tests to determine significance (α=0.10) are shown above the plotted 

points. Information in brown applies to the 14-year trends and black applies to 20-year trends.

Figure 25. Trend slopes by geologic class for the seven metrics with longterm trend slopes indicating significant geologic effects. The p-values from 

one-sample t-tests to determine significance (α=0.10) are shown above the 

plotted points. Information in brown applies to the 14-year trends and black 

applies to 20-year trends.
44 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Likewise, for the seven AM metrics whose trend slopes 

depended on geology, the 20-year timeframe was used for the 

same reason; that is, maximum observed trend slopes among 

sites tended to be highest for the 20-year timeframe. Thus, 

these results report the “maximum potential” power to detect 

change rather than the average. A type I effort rate of 0.10 and 

a power of 0.20 were used for all calculations. As with trend 

analyses, power calculations were made on transformed data 

where appropriate (table 10). Direct comparisons between the 

two groups should be made with caution because variances 

used in MDC calculations were computed using different 

sample sizes. Specifically, sample sizes for the metrics whose 

trends depended on geology were smaller than those used 

for park-wide metrics because N equaled the number of sites 

within a geologic type rather than total number of sites. Thus, 

one would expect lower power (and therefore greater YTDC) 

for those metrics. In addition, calculations were made for a 

given metric only for geologic types that showed statistically 

significant trend slopes (table 13).

In terms of metrics having park-wide trends, those for 

total richness and EPT-r were strongest. Given the average 

10-year trend slope, park-wide changes in these metrics could 

be detected in as few as three years (fig. 27). The Leuctra:P% 

metric trend also was strong, with significant differences 

detectable in as few as six years. The DOM2 and SimpD metric trends were weaker and required between 15 and 18 years 

of data to detect a significant change.

For those metrics whose trend varied with geology, at 

least 20 years of data were typically required to detect a significant change (fig. 27), indicating the weakness of observed 

trends and highlighting the importance of long-term monitoring data. Additionally, those metrics whose trends were indicative of declining stream condition, which always occurred in 

the siliciclastic watersheds, could be detected in fewer years 

than corresponding metrics indicative of improving stream 

condition, which always occurred in granitic and basaltic 

watersheds (fig. 27). 

0 1,000 2,000 3,000 –0.10

–0.05

0

0.05

0.10

0 1,000 2,000 3,000

–0.06

–0.04

–0.02

0

0.02

0 1,000 2,000 3,000 –0.6

–0.3

0

0.3

0.6

0 1,000 2,000 3,000

–0.06

–0.04

–0.02

0

0.02

Watershed area, in hectares Watershed area, in hectares

R2 =0.174 R2 =0.268

R2

 =0.141

R2

 =0.204

Figure 26. Graphs showing trend slope as a function watershed area for 

selected aquatic macroinvertebrate metrics.

EPT%

PTV Intol%

Shredder%

Trend slopes

Figure 26. Trend slope as a function watershed area for selected aquatic 

macroinvertebrate metrics. 
Analysis of Aquatic Macroinvertebrate Data 45

Siliciclastic

Geologic class

EXPLANATION

Granitic

Basaltic

*

*

*

*

*

*

A. Metrics showing park-wide trends

Richness EPT-r LeuctraP

Aquatic macroinvertebrate metric

Aquatic macroinvertebrate metric

DOM2 SimpD

Number of years required to detect change (YTDC)

0

4

8

12

16

20

B. Metrics with trends that depend on geology

EPT% HydroT PTV Intol Scrap Shred SCI

0

10

20

30

40

50

60

Figure 27. Graphs showing the number of years required to detect a 

statistically significant change in metric means for A, metrics that 

showed significant trends park-wide, and B, metrics with trends 

dependent on geology. Asterisks above the bars indicates that trends 

were positive (improvement in stream condition). All other trends 

were negative (declining stream condition).

Figure 27. The number of years required to detect a statistically significant 

change in metric means for A, metrics that showed significant trends parkwide, and B, metrics with trends dependent on geology. Asterisks above the 

bar indicates that trends were positive (improvement in stream condition). 

All other trends were negative (declining stream condition). 
46 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Analysis of Fish Data

Stream fishes are a valued biological resource within the 

park, and a robust fish-community dataset has been collected 

by NPS personnel over several years. This dataset was used 

to evaluate the status and trends of stream fish communities 

within the study area. Analyses focused on geology and watershed size as environmental predictors of stream fish communities because of their demonstrated importance in relation to 

acidic deposition (Bulger and others, 1995) and fish community composition (Roberts and Hitt, 2010).

Fish Data Structure and Analysis

Park personnel have collected stream fish data since 

1982, with an emphasis on understanding brook-trout population dynamics (Atkinson, 2002). The 1996–2009 fish dataset 

includes 107 sample sites and 512 individual survey events 

(fig. 28; table 14). This dataset provides an unparalleled resolution for spatial and temporal analysis of Appalachian stream 

fish communities.

Stream sites were sampled using 3-pass backpack 

electrofishing techniques within standardized 100-m reaches. 

Captured fishes were identified to species, counted by pass 

number, and released. The majority of sampling events 

occurred during June and July, but some samples were collected during May, August, and September (table 14). The 

number of sampling sites per year has generally increased over 

time, as has the variation in sample size (table 14). Because 

the proportion of samples per month remained nearly constant 

over time (table 14), however, temporal changes are unlikely 

to be confounded by sampling effort.

Counts of brook trout young-of-year (YOY) were distinguished from adults (age greater than 1 year; age 1+) (Atkinson, 2002). Although it is possible that some juvenile fishes 

were erroneously considered “adults” in this analysis, brook 

trout spawning may occur in either the first or second full 

year in the study area (Jenkins and Burkhead, 1994), so it was 

appropriate to distinguish YOY from adult (age 1+) fish in this 

analysis. The total abundances were estimated for each species 

as the sum of abundances within each pass.

Subsets of the fish dataset were evaluated to accomplish 

different objectives. First, all sites were evaluated to describe 

the approximate relative abundance of all fish species in the 

study area (n=107). Second, sites that supported at least three 

species and over 9 years of annual survey data (n=25) were 

selected for community trend analyses. Third, all sites that 

supported more than 9 years of annual survey data (n=32) 

were selected for brook-trout trend analyses; no required species richness threshold was necessary in the selection of these 

sites. Of the 32 brook-trout focal sites, 22 (69 percent) also 

were included as fish community sites. Community and brook 

trout focal sites included locations throughout the study area 

(fig. 28).

Fish-community and brook-trout population analyses 

followed similar methods for water quality and AM communities. Specifically, simple-linear-regression slope coefficients 

were calculated for biotic metrics (fish richness, brook trout 

YOY, and adult brook-trout abundances) as a function of year. 

In addition, the per-capita intrinsic rate of increase (r) was 

calculated as 

r=ln(Nt/Nt-1) (1)

where

 N is the adult abundance, and

t is the sampling year.

This approach provided estimates of increasing and 

decreasing trends for fish metrics over time. Linear models for 

temporal trend analysis were evaluated using the same time 

periods, significance criteria, and geologic classes used for the 

water-quality and AM trend analyses. GLM was used to evaluate the importance of geologic class, upstream watershed area, 

and interactive effects on fish metrics (species richness, YOY, 

and adult brook-trout abundance) among time periods.

Table 14. Number of fish community sampling events per month 

from 1996 to 2009.

Year

Month

May June July August September Total

1996 0 10 13 6 0 29

1997 0 14 14 4 0 32

1998 0 16 18 0 0 34

1999 0 17 19 0 0 36

2000 0 18 17 2 0 37

2001 0 17 19 4 0 40

2002 0 20 20 1 0 41

2003 0 16 20 0 0 36

2004 0 21 18 10 0 49

2005 0 8 7 0 0 15

2006 3 17 25 16 0 61

2007 0 2 4 0 0 6

2008 1 26 29 7 0 63

2009 0 12 14 5 2 33

Total 4 214 237 55 2 512
Analysis of Fish Data 47

!

!

!

!

!!

!

!

!

!

!

!

!

!

!

!

!

!

!!

!

! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!! ! ! !

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!

!!

!

!

!

!

! !

!

!

!

!

!

!

!

!

!

! !

!!!

!

!

!!

!

!

!

!

!

!

!

!

!

!

!

!

!

! !

!

!

!

!

!!

!

!

!

!!

!!

!

! !

!

!

!!!

!

!

!

! !

!

!

!

!

!!

!

!

!

!!

!!

!

! !

! !

!

!!!

!

!

! !

!!

! !

!

!

10 20 KILOMETERS

5 10 20 MILES

EXPLANATION

All fish sites (1996–2009)

Brook trout and community focal site

Fish sampling sites

Community focal site only

Brook trout focal site only

N 0

0

Figure 28. Map showing all fish sampling site locations, and those selected for brook trout analyses and 

community analyses, within Shenandoah National Park. Location of study area provided in figure 1.

Shenandoah National Park 

 boundary

Streams

Siliciclastic

Bedrock geology

Granitic

Basaltic

Figure 28. All fish sampling site locations, and those selected for brook trout analyses and community 

analyses, within Shenandoah National Park. Location of study area provided in figure 1.
48 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Fish Community Characteristics

The park’s fish dataset contained records of 261,028 

individual fish and 33 species (n=107). Members of the family 

cyprinidae were numerically dominant (table 15), primarily 

because of the large number of blacknose dace (Rhinichthys 

atratulus) in the dataset (37 percent, table 15). Salmonid fishes 

were second-most abundant in the dataset, primarily because 

of brook trout abundances (25 percent), rather than introduced 

rainbow trout (Oncorhynchus mykiss) and (or) brown trout 

(Salmo trutta) (table 15). Cottid and catostomid fishes collectively contributed approximately one-fifth of the total individuals in the dataset because of relatively high abundances of 

mottled sculpin (Cottus bairdi) and torrent sucker (Thoburnia 

rhothoeca) (table 15). Most of the taxonomic diversity was represented by relatively low abundances. For example, 22 species 

(67 percent) were represented by fewer than 1,000 individuals, 

and 12 species (36 percent) contained fewer than 100 individuals in the dataset. In contrast, the most abundant species (blacknose dace) consisted of nearly 100,000 individuals (table 15).

Table 15. Fish species and abundances reported from 1996 to 2009 from 3-pass backpack electrofishing surveys at 107 sites. 

[Species codes correspond to those shown in figure 29. %, percentage; <, less than]

Family (% of total) Scientific name Common name Code Number of 

individuals % of total

Anguillidae (0.2) Anguilla rostrata American eel AME 428 0.16

Catostomidae (9.4)

Catastomus commersoni White sucker WHS 496 0.19

Hypentelium nigricans Northern hog sucker NHS 85 0.03

Thoburnia rhothoeca Torrent sucker TOS 24,033 9.21

Centrarchidae (0.1)

Ambloplites rupestris Rock bass ROB 30 0.01

Lepomis auritus Redbreast sunfish RDB 8 <0.01

Lepomis macrochirus Bluegill BLG 2 <0.01

Lepomis gibbosus Pumpkinseed PKS 3 <0.01

Micropterus dolomieu Smallmouth bass SMB 322 0.12

Micropterus salmoides Largemouth bass LMB 2 <0.01

Cottidae (11.3)

Cottus bairdi Mottled sculpin MTS 28,705 11.00

Cottus girardi Potomac sculpin POS 653 0.25

Cyprinidae (49.0)

Campostoma anomalum Central stoneroller CES 3,591 1.38

Clinostomus funduloides Rosyside dace RYD 8,740 3.35

Cyprinella analostana Satinfin shiner SNS 3 <0.01

Exoglossum maxillingua Cutlips minnow CUM 133 0.05

Luxilus cornutus Common shiner COS 1,226 0.47

Nocomis leptocephalus1 Bluehead chub BLC 4,904 1.88

Phoxinus oreas Mountain redbelly dace MRD 2,335 0.89

Pimephales notatus Bluntnose minnow BLM 36 0.01

Rhinichthys atratulus Blacknose dace BKD 95,624 36.63

Rhinichthys cataractae Longnose dace LGD 10,829 4.15

Semotilus atromaculatus Creek chub CRC 255 0.10

Semotilus corporalis Fallfish FAF 254 0.10

Ictaluridae (0.1)

Ameiurus natalis Yellow bullhead YEB 6 <0.01

Ameiurus nebulosus Brown bullhead BRB 1 <0.01

Noturus insignis Margined madtom MAM 296 0.11

Percidae (4.7)

Etheostoma blennioides Greenside darter GRD 3 <0.01

Etheostoma flabellare Fantail darter FAD 11,945 4.58

Etheostoma olmstedi Tessellated darter TED 414 0.16

Salmonidae (25.2) Oncorhynchus mykiss Rainbow trout RBT 56 0.02

Salmo trutta Brown trout BRT 367 0.14

Salvelinus fontinalis Brook trout BKT 65,240 24.99

1

Combined with records of river chub (Nocomis micropogon) due to taxonomic uncertainty for juvenile specimens.
Analysis of Fish Data 49

Temporal Trends and Spatial Patterns in 

Fish Communities

Most species did not exhibit dramatic changes in average 

abundance during 1996–2009, such as park-wide colonization 

or extirpation (fig. 29). Some patterns, however, merit further 

investigation: potential increases in average abundance of 

American eel (Anguilla rostrata; AME), margined madtom 

(Noturus insignis; MAM), mountain redbelly dace (Phoxinus oreas; MRD), and decreases in average abundances of 

creek chub (Semotilus atromaculatus; CRC) and the common 

shiner (Luxilus cornutus; COS) (fig. 29), although variation in average abundance may have been influenced by the 

number of sampling events per year (table 14). For example, 

low abundances during 2007 were probably a reflection of the 

relatively small number of samples that year, and increases in 

American eel abundance may be attributed to the removal of 

Embrey Dam on the Rappahannock River in 2004 (Hitt and 

others, 2012).

Fish species richness in focal sites tended to increase or 

remain constant over time; no sites exhibited an obvious loss 

of richness (fig. 30). Temporal increases in richness were most 

pronounced in the North Fork of Thornton River (1F030) and 

in Hughes River (2F038), where richness appeared to increase 

from approximately 9 to 13 species and from 11 to 15 species, respectively (fig. 30). Richness at sites supporting three 

species in 1996 (typically brook trout, blacknose dace, and 

sculpin) remained relatively constant over time. In contrast, 

sites with four or more species in 1996 showed more variation 

in richness over time (sites 2F039, 2F040, 2FVA4). Collectively, these results suggest that fish dispersal and recolonization dynamics may be operating across larger spatial scales 

than the 100-meter sampling reach (Brown and others, 2011).

Mean species richness was greater in granitic and basaltic 

watersheds than in siliciclastic watersheds and increased with 

watershed area (fig. 31). Granitic watersheds supported sites 

having the greatest fish richness, with the opposite being true 

for siliciclastic watersheds (fig. 31). GLM results indicated 

watershed area and geologic class affected mean species richness within all time periods (5, 10, and 14 years) (table 16). 

No interactive effects of geology and watershed area on mean 

species richness were observed (table 16).

Species richness tended to increase over time, but these 

changes were evident only for the 14-year analysis period 

(fig. 32). In the 14-year analysis, geologic classes exhibited 

important differences: sites draining granitic and basaltic 

watersheds showed increasing species richness, whereas 

sites draining siliciclastic watersheds showed no such trends 

(fig. 32). In addition, larger streams tended to accumulate 

more species over time than smaller streams. Although the 

effect of watershed area on fish richness trends was largely 

insignificant in the GLMs (table 16), this was primarily 

because of decreasing species richness observed at one site 

Table 16. Stepwise general linear model results for effects of geology and watershed area on stream fish community richness, 

brook trout young-of-year abundance, and brook trout adult (age greater than 1 year: age 1+) abundance among time periods of 5,

10, and 14 years. Non-significant terms were sequentially removed from models. 

[ns, not significant for α=0.10; YOY, young-of-year]

Response Metric Years Geologic class Watershed area Interaction

F-statistic p-value F-statistic p-value F-statistic p-value

Mean

Richness 5 3.72 0.043 14.50 <0.001 ns ns

10 3.60 0.047 16.99 <0.001 ns ns

14 3.83 0.040 18.39 <0.001 ns ns

YOY 5 5.76 0.009 ns ns ns ns

10 7.77 0.003 4.90 0.036 ns ns

14 8.49 0.001 4.86 0.037 ns ns

Adult 5 ns ns 7.47 0.011 ns ns

10 2.94 0.071 5.94 0.022 ns ns

14 4.55 0.020 7.622 0.010 ns ns

Slope

Richness 5 ns ns ns ns 3.09 0.069

10 ns ns ns ns ns ns

14 ns ns ns ns ns ns

YOY 5 4.27 0.022 ns ns ns ns

10 ns ns ns ns ns ns

14 ns ns ns ns ns ns

Adult 5 6.47 0.005 ns ns ns ns

10 ns ns ns ns ns ns

14 3.93 0.032 ns ns ns ns
50 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

● ● ● ● ● ● ● ● ● ● ●

●

● ● ●

●

●

●

● ● ● ● ● ●

● ● ● ● ● ● ●

●

● ● ● ●

●

●

● ● ● ● ●

● ● ● ● ●

● ●

● ● ●

●

●

● ●

●

● ●

● ● ● ● ●

●

● ● ● ● ● ●

● ● ● ● ● ● ●

●

● ● ● ●

● ● ● ● ● ●

●

●

● ● ● ●

●

●

● ● ● ● ● ● ● ● ● ●

● ●

● ● ● ●

●

●

●

● ●

● ● ● ● ● ● ● ● ●

●

●

● ●

●

● ●

●

●

● ●

● ●

● ●

● ● ● ●

● ●

●

●

●

●

● ● ● ● ● ●

●

●

●

●

●

●

● ●

● ● ●

●

● ● ● ● ● ● ● ● ● ●

●

●

●

● ● ●

● ●

●

● ● ● ● ● ● ●

● ●

● ● ● ● ● ● ● ● ● ●

●

●

● ●

● ●

● ● ● ● ● ● ● ● ●

● ● ●

●

●

● ●

●

● ● ● ●

●

● ●

●

●

● ● ● ●

●

●

● ● ● ● ● ● ● ● ●

●

● ●

● ●

● ●

● ● ● ● ●

●

● ● ●

●

●

●

● ●

● ● ●

● ● ●

●

●

● ● ●

●

● ● ● ● ●

Average species abundanceAverage species abundance

Figure 29. Graphs showing the temporal variation in average abundance of species within Shenandoah National 

Park community focal sites (n=25). 

0

2

4

6

8

10

12 AME

0

20

40

60

80 BHC

0

200

400

600

800 BKD

0

100

200

300

400

500 BKT

0

2

4

6

8

10

12

BRT

0

100

200

300 CES

0

5

10

15

20

25

30 COS

0

1

2

3

4

5

CRC

0

1

2

3

4

5

6

7 CUM

0

20

40

60

80

100

120

140

FAD

0

5

10

15

20

25

FAF

0

50

100

150

200

250 LGD

0

5

10

15

20

25

30

MAM

0

10

20

30

40

50

60

70

80

MRD

0

50

100

150

200

250

300

350 MTS

0

1

2

3

4

5

RBT

0

50

100

150

200

RYD

0

2

4

6

8

10

12

14

SMB

0

2

4

6

8

10

TED

0

200

400

600

800

1996 1998 2000 2002 2004 2006 2008 2010

TOS

0

1

2

3

4

1996 1998 2000 2002

Year

Year

2004 2006 2008 2010

WHS

EXPLANATION

1 standard deviation above mean

● ● Mean

Figure 29. The temporal variation in average abundance of species within Shenandoah National Park community 

focal sites (n=25). 
Analysis of Fish Data 51

●

●

●

●

● ● ● ● ●

● ●

● ● ● ● ● ● ● ● ● ● ●

● ● ●

●

●

●

●

●

● ●

● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ●

●

●

●

●

●

● ● ● ● ● ● ● ● ● ● ●

● ●

●

● ●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

● ● ●

●

●

●

●

● ● ● ●

● ● ● ● ● ● ● ● ● ● ●

●

● ● ●

●

● ● ● ● ● ● ● ●

●

● ●

●

● ● ●

●

●

●

● ●

● ●

●

●

● ●

●

●

●

● ● ● ● ●

● ● ● ●

●

●

●

●

● ●

●

●

●

●

●

●

●

●

●

●

●

● ●

● ●

●

●

●

●

●

●

● ● ● ● ●

●

●

● ● ● ● ●

● ●

● ●

● ●

●

●

● ●

●

●

● ● ● ● ● ● ● ● ● ● ● ● ●

● ●

● ● ●

●

● ●

● ● ●

●

●

●

●

●

● ● ●

●

● ● ● ●

●

● ● ● ● ● ● ● ●

●

● ● ● ●

● ● ● ● ● ● ● ● ●

●

● ● ●

● ● ● ● ● ● ● ● ● ● ● ● ●

● ● ● ● ● ● ● ● ● ● ●

Figure 30. Graphs showing the temporal variation in fish species richness for community focal sites by National 

Park Service (NPS) site code. Site locations and attributes are provided in table 9. Fish species richness

Fish species richness

Year

Year

4

5

6

7

8

3

4

5

6

7

9

10

11

12

13

14

2.5

3.5

4.5

5.5

2.0

2.5

3.0

3.5

4.0

4.5

2.5

3.5

4.5

5.5

5.0

5.5

6.0

6.5

7.0

7.5

3.0

3.5

4.0

4.5

5.0

5.5

6.0

6.5

5.0

5.5

6.0

6.5

7.0

7.5

1.5

1.0

2.0

2.5

11

12

13

14

15

16

4.0

4.2

4.4

4.6

4.8

5.0

5.2

3.0

3.5

4.0

4.5

5.0

5.5

2.0

2.5

3.0

3.5

4.0

4.5

10

12

14

16

18

20

6

7

8

9

10

12

4.0

4.5

5.0

5.5

6.0

6.5

3.0

3.2

3.4

3.6

3.8

4.0

4.2

3

4

5

6

7

4.0

4.5

5.0

5.5

6.0

6.5

7.0

7.5

9.0

9.5

10.0

10.5

11.0

11.5

12.0

12.5

3.0

3.2

3.4

3.6

3.8

4.0

4.2

1.0

1.5

2.0

2.5

3.0

3.0

2.0

2.5

3.0

3.5

4.0

2.0

2.5

3.0

3.5

4.0

1F003 1F007

1F030 1F118

1F132

1FVA1

1FVA2

1FVA3

2F016

2F017

2F038

2F039

2F040 2F072

2F093 2F095

2F135

2FVA4

3F044 3F045

3F084 3F123

3F124

1996 1998 2000 2002 2004 2006 2008 2010

3F125

1996 1998 2000 2002 2004 2006 2008 2010

3F126

Figure 30. The temporal variation in fish species richness for community focal sites by National Park Service (NPS) site 

code. Site locations and attributes are provided in table 9.
52 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

●

●

●

●●

●

●

●

●

●

●

●

●

●

●

●

●

● ● ●

●

●

●

●

●

●

●

● ● ●●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ● ●●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

● ●

●

● ● ●●

EXPLANATION

90th percentile

75th percentile

Median

25th percentile

10th percentile

Outlier

5−year mean

0

0

0

2

4

6

8

10

12

14

16

16 10−year mean

2

4

6

8

10

12

14

14−year mean

2

4

6

8

10

12

14

Siliciclastic Granitic Basaltic

Fish richness

0

0

0

2

4

6

8

10

12

14

16

16

2

4

6

8

10

12

14

2

4

6

8

10

12

14

Fish richness

Geologic class

5−year mean

10−year mean

14−year mean

500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

Watershed area, in hectares

Figure 31. Graphs showing distributions of mean fish richness by geologic class and against watershed 

area for three time periods analyzed. Letters on boxplots indicate significant differences between 

classes (α = 0.10). A least-squares regression line is shown for watershed-area graphs.

A

A, B

B

A B A, B

A B A, B

Figure 31. Distributions of mean fish richness by geologic class and against watershed area for 

three time periods analyzed. Letters on boxplots indicate significant differences between classes 

(alpha = 0.10). A least-squares regression line is shown for watershed-area graphs. 
Analysis of Fish Data 53

●

●

●

●

●

●

●

●

●

●

●

●

●

●

● ●

●

●

●

●

●

●

●

EXPLANATION

●

Not significant (α = 0.10)

Trend slope

Significant (α = 0.10)

5−year slope

10−year slope

14−year slope

5−year slope

10−year slope

14−year slope

–0.4

–0.2

0

0.2

0.4

0.6

0.6

0.6

–0.4

–0.2

0

0.2

0.4

–0.4

–0.2

0

0.2

0.4

Siliciclastic Granitic Basaltic

Slope of fish richness over time

−0.4

−0.2

0.0

0.2

0.4

0.6

0.6

0.6

−0.4

−0.2

0.0

0.2

0.4

−0.4

−0.2

0.0

0.2

0.4

Slope of fish richness over time

Geologic class

500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

Watershed area, in hectares

Figure 32. Graphs showing fish richness trend slope as functions of geologic class and watershed 

area for the three time periods analyzed. 

Figure 32. Fish richness trend slope as functions of geologic class and watershed area for the 

three time periods analyzed. 
54 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

(2F093). Removal of this site from the analysis yielded a 

significant relation between watershed area and species richness slopes (p<0.01). The outlier site 2F093 is notable for its 

history: it was established in 1996 after a major flood in 1995, 

which reworked the substrate and removed riparian canopy 

cover. The decline in species richness for this site suggests 

that floods may have long-term effects on fish community 

composition, despite the recovery of brook trout (Roghair and 

Dolloff, 2005) and AMs (Snyder and Johnson, 2006). Another 

fish community focal site, however, was scoured by the 1995 

flood (2F072; Staunton River), and this site showed increasing 

species richness (14-year slope coefficient=0.172; watershed 

area=1,065 ha). Some increases in richness were associated 

with the appearance of native American eel at a site, whereas 

other increases were associated with the occurrence of nonnative rainbow trout or potentially non-native torrent sucker.

Detection of temporal trends in fish richness required 

more than 5 years of data. No significant (p<0.10) trend slopes 

were detected for the 5-year period, whereas 6 significant 

trends were detected for the14-year period (fig. 32). Statistical significance in this regard generally was detected only 

for positive or negative changes in slope greater than 0.10, 

suggesting the sensitivity of linear regression models to detect 

a net gain or loss of 1 species every 10 years, on average. This 

finding indicates the necessity of long-term monitoring data 

for detecting changes in stream-fish richness.

Temporal Trends and Spatial Patterns in Brook 

Trout Population

Mean brook-trout YOY abundances were substantially 

greater in granitic and basaltic watersheds than in siliciclastic 

watersheds, but were invariant with respect to watershed size 

(table 16; fig. 33). The effect of geologic class on mean YOY 

abundance was observed in the 5-year analysis, as well as 

the 10- and 14-year analyses (table 16; fig. 33). Siliciclastic 

watersheds supported approximately half of the mean YOY 

abundance observed in other geologic classes (fig. 33).

Temporal trends in YOY abundance varied substantially 

among time periods and geologic classes (fig. 34). The 5-year 

time periods generally showed more variation in YOY abundance trends than longer time intervals. For the 10- and 14-year 

analyses, more sites were characterized by decreasing trends in 

YOY abundance than increasing trends (fig. 34). Basaltic sites 

generally showed the fastest rate of decrease (largest negative 

slope coefficients), but these declines mostly were not significant (fig. 34). Stream volume was probably not important for 

YOY abundance trends, because decreasing YOY abundances 

were observed across sites having small and large upstream 

watershed areas (fig. 34). The apparent increases in YOY 

abundance within granitic sites for the 5-year analysis were 

not observed in analyses of longer timespans (fig. 34).

As with YOYs, the mean abundance of adult (age 1+) 

brook trout was greater in granitic and basaltic watersheds 

than in siliciclastic watersheds (table 16; fig. 35). Larger 

watersheds tended to support lower average abundances of 

adult brook trout (fig. 35), and these differences were statistically significant (table 16). Five-year temporal trends in adult 

abundances were more variable than trends for longer time 

spans (fig. 36). Unlike YOY trends, however, trends in adult 

abundances at some sites increased over time in the 14-year 

analysis. These increasing trends were more prevalent for 

siliciclastic and granitic sites than basaltic sites, where both 

increasing and decreasing trends were observed in relatively 

equal amounts (fig. 36). Watershed size was not associated 

exclusively with increasing or decreasing trends in adult 

abundance over time (fig. 36). Model sensitivity in detecting 

significant temporal change in adult brook-trout abundances 

increased with the length of the period analyzed, underscoring 

the importance of long-term monitoring data.

Per capita intrinsic rates of increase (r) for adult brook 

trout populations generally centered around zero, indicating population stability over the period of record (fig. 37). 

Preliminary analyses revealed insignificant differences among 

geologic classes or watershed size, but other patterns merit 

additional consideration. Several sites exhibited decreasing variation in r over time. For example, intrinsic rates of 

increase in Jeremy’s Run (1F007) and Hughes River (2F038) 

showed more variability in the first half of the time series than 

the second (1997–2003 versus 2003–09) (fig. 37). This result 

suggests possible differences in the natural disturbance regime 

among park locations. In addition, all focal sites showed 

consistent patterns of change during selected periods. For 

example, of the 22 sites with data, 19 (86 percent) showed 

r increasing from 1997 to 1998 and decreasing the following year. This result suggests an important role of densityindependent environmental factors, such as flow regime, apart 

from the density-dependent effects of predation and competition. Further research is needed to understand whether or not 

the temporal synchrony observed here is also observed across 

larger regions of the southern Appalachians.

This analysis provided several inferences regarding fish 

community and brook-trout population dynamics within the 

park. Fish species richness was found to be depressed in acidsensitive (siliciclastic) watersheds, as reported by Bulger and 

others (1995), and fish species richness has increased over time

in granitic and basaltic watersheds not sensitive to acidification. Likewise, abundances of YOY and adult brook trout were 

depressed in siliciclastic watersheds, but in contrast to fish 

species richness, increasing trends in adult brook trout abundance were detected in acid-sensitive sites over the 14-year 

analysis timeframe. The observed long-term increases in adult 

brook trout abundance in acid-sensitive streams do not reflect 

improvements in fish community composition (species richness). Analysis of brook-trout population growth rates revealed 

stable populations (r ~ 0) and highly synchronous patterns 

of annual variation park-wide. Further analysis of spatially 

structured environmental factors, such as streamflow and water 

temperature, could therefore enable predictions of brook trout 

responses to climate change within the park and throughout the 

Appalachian Mountains. Conservation of fish biodiversity and 

brook trout populations in acid-sensitive watersheds remains an 

important challenge for research and management.
Analysis of Fish Data 55

EXPLANATION

90th percentile

75th percentile

Median

25th percentile

10th percentile

Outlier

500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

Watershed area, in hectares

Granitic

5−year mean

0

50

100

150

10−year mean

0

50

100

150

14−year mean

0

50

100

150

200

200

200

Siliciclastic Basaltic

Geologic class

Brook trout YOY abundance

0

50

100

150

0

50

100

150

0

0

50

100

150

200

200

200

Brook trout YOY abundance

5−year mean

10−year mean

14−year mean

Figure 33. Graphs showing distributions of brook trout young-of-year (YOY) mean abundance by 

geologic class and against watershed area for three time periods analyzed. Letters on boxplots 

indicate significant differences between classes (α = 0.10). Watershed area plots show a 

least-squares regression line.

A A, B

B

A B B

A B B

Figure 33. Distributions of brook trout young-of-year (YOY) mean abundance by geologic 

class and against watershed area for three time periods analyzed. Letters on boxplots indicate 

significant differences between classes (alpha = 0.10). Watershed area plots show a least-squares 

regression line. 
56 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

●

●

●

●

●

● ●

●

●

EXPLANATION

Not significant (α = 0.10)

Trend slope

Significant (α = 0.10)

Watershed area, in hectares

0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

Figure 34. Graphs showing brook trout young-of-year (YOY) abundance trend slopes by geologic 

class and against watershed area for three time periods analyzed.

5−year slope

–50

–100

0

50

100

–100

100

–100

100

10−year slope

–50

0

50

14−year slope

–50

0

50

Siliciclastic Granitic Basaltic

Slope of fish richness over time

–50

–100

0

50

100

–100

100

–100

100

–50

0

50

–50

0

50

Slope of fish richness over time

Geologic class

5−year slope

10−year slope

14−year slope

Figure 34. Brook trout young-of-year (YOY) abundance trend slopes by geologic class and against 

watershed area for three time periods analyzed. 
Analysis of Fish Data 57

EXPLANATION

90th percentile

75th percentile

Median

25th percentile

10th percentile

Outlier

0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

5−year mean

0

50

100

150

200

200

200

10−year mean

0

50

100

150

14−year mean

0

50

100

150

Siliciclastic Granitic Basaltic

Geologic class

Brook trout adult (age I+) abundance

0

50

100

150

200

200

200

0

50

100

150

0

50

100

150

Brook trout adult (age I+) abundance

5−year mean

10−year mean

14−year mean

Watershed area, in hectares

Figure 35. Graphs showing distributions of brook trout adult (age greater than 1 year; age 1+) mean 

abundance by geologic class and against watershed area for three time periods analyzed. Letters 

on boxplots indicate significant differences between classes (α = 0.10). Watershed area plots 

show a least-squares regression line.

A A, B B

A A A

A A A

Figure 35. Distributions of brook trout adult (age greater than 1 year; age 1+) mean abundance by 

geologic class and against watershed area for three time periods analyzed. Letters on boxplots indicate 

significant differences between classes (alpha = 0.10). Watershed area plots show a least-squares 

regression line. 
58 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

EXPLANATION

Not significant (α = 0.10)

Trend slope

Significant (α = 0.10)

Figure 36. Graphs showing brook trout adult (age greater than 1 year; age 1+) abundance trend 

slopes by geologic class and against watershed area for three time periods analyzed. 

5−year slope

–20

–30

–40

–10

0

10

20

30

40

–30

–40

30

40

–30

–40

30

40

10−year slope

–20

–10

0

10

20

14−year slope

–20

−10

0

10

20

Siliciclastic Granitic Basaltic

Slope of brook trout adult (age I+) abundance over time

–20

–30

–40

–10

0

10

20

30

40

–30

–40

30

40

–30

–40

30

40

–20

–10

0

10

20

–20

–10

0

10

20

Slope of brook trout adult (age I+) abundance over time

Geologic class

5−year slope

10−year slope

14−year slope

500 1,000 1,500 2,000 2,500 3,000 3,500 4,000

Watershed area, in hectares

0

Figure 36. Brook trout adult (age greater than 1 year; age 1+) abundance trend slopes by 

geologic class and against watershed area for three time periods analyzed. 
Analysis of Fish Data 59

–2

–1

0

1

2

3

3

1F003 1F007 1F030

1F118

–2

–1

0

1

2 1F132 1F138

–2

–1

0

1

2

1FVA1 1FVA2 1FVA3

–2

–1

0

1

2

3

4

3

3

3

3

3

4

3

4

4

2F016 2F017 2F038

–2

–1

0

1

2

2F039 2F040 2F072

–2

–3

–1

0

1

2 2F074 2F075 2F093

–2

–1

0

1

2

2F095 2F135

2FVA4

–2

–1

0

1

2 3F044 3F045 3F102

–2

–1

0

1

2 3F103 3F107 3F109

–2

–1

0

1

2

3F123 3F124

1996 2000 2004 2008

3F125

–2

–1

0

1

2

–2

–1

0

1

2

3

3

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

3

4

3

3

3

3

3

4

3

4

4

–2

–1

0

1

2

–2

–3

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

3

3

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

3

4

3

3

3

3

3

4

4

–2

–1

0

1

2

–2

–3

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

–2

–1

0

1

2

1996 1998 2000 2002 2004 2006 2008 2010 1998 2002 2006 2010

1998 2002 2006 2010

3F126

1996 2000 2004 2008

3F128

Per-capita intrinsic rate of increase (r)

Figure 37. Graphs showing per capita rate of increase (r) for adult brook trout at 32 focal sites.

Year Year

Year

Figure 37. Per capita rate of increase (r) for adult brook trout at 32 focal sites.
60 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Relations Among Environmental 

Factors, Water Quality, Aquatic 

Macroinvertebrates, and Fish 

Similar temporal trends among water-quality, AM, and 

fish metrics were evaluated in order to derive hypotheses 

about causal explanations for temporal trends. Unfortunately, 

the relations cannot be modeled directly because data were not 

collected consistently from common sites. Three anthropogenic stressors are most likely to affect stream condition in the 

park: acidic deposition, pest-induced declines in forest health 

(primarily caused by the gypsy moth), and climate change. 

Local factors such as human development are less important 

in the park because of the protection afforded by its status as a 

national park unit. 

Stream Acidification

Federal mandates associated with Title IV amendments 

of the Clean Air Act established in 1990 have significantly 

improved air quality nationally, including substantial reductions in sulfur dioxide (SO2

) and nitrous oxide (NOx) emissions, the primary gasses that contribute to acidic deposition. In 

a recent national acidic deposition assessment, Burns and others (2011) found that in addition to substantial improvements 

in air quality, there is evidence that the chemistry of some 

lakes and streams are recovering, although these improvements 

have been modest and sparse. In particular, streams in the 

mid-Atlantic have shown minimal evidence of improvement, 

with only 12 percent of monitored streams showing increasing 

trends in ANC and decreasing trends in SO4

2– concentrations 

between 1990 and 2008 (Burns and others, 2011). 

The assessment of stream-water quality in the park 

presented herein shows a conflicting pattern of acidity trends. 

Long-term trends in ANC and pH suggest some improvement 

in terms of decreasing acidity, but short-term trends suggest 

increasing acidity. Long-term increases in pH occurred across 

all geology classes park-wide. In contrast, the long-term 

increases in ANC were largely limited to sites in watersheds 

underlain by basaltic geology (figs. 7 and 8). Decreasing trends

in SO4

2– concentrations indicated long-term improvements in 

basaltic watersheds, but increasing trends in siliciclastic and 

granitic watersheds indicate deteriorating conditions in those 

watersheds (fig. 9). These patterns conform to the general 

conclusion of Burns and others (2011) that recovery of midAtlantic streams to acidic deposition is weak at best. 

There is some evidence that the AM community may 

have responded favorably to the long-term improvements 

in stream pH observed park-wide. Specifically, the proportion of the acid tolerant stonefly Leuctra (Leuctra:P% metric) 

declined park-wide, a result consistent with park-wide declines 

in stream acidity. The Leuctra:P% metric was selected for 

use in the long-term biological monitoring program in the 

hope that it would be a stress-specific indicator of acidity in 

streams;, however, Leuctra abundance is also strongly related 

to water temperature. In contrast to the Leuctra:P% metric, 

however, the majority of AM and water-quality data suggest water quality is not improving, at least not enough to 

positively affect AMs. First, the trend slopes for long-term 

increases in pH were modest (mean=approximately 0.005 per 

year, fig. 8). Thus, on the average, pH increased about 0.1 pH 

unit over 20 years, an improvement not likely to be biologically significant. Second, short-term trends in pH actually 

showed the opposite trend (decreasing pH; fig. 8), whereas 

the Leuctra:P% metric declined within both long-term and 

short-term timeframes (fig. 24). Finally, other AM metrics that 

would also be expected to improve in response to decreasing 

acidity or increased water quality in general (EPT%, PTV, 

Intol%, Scraper%, and Hydro:T%), instead, were indicative 

of long-term stream degradation (fig. 25), and in many cases 

(EPT%, PTV, Scraper%, Hydro:T%, and Intol%; fig. 25) 

the evidence for degradation was most pronounced in those 

streams draining siliciclastic bedrock. Thus, in contrast to 

long-term trends in the Leuctra:P% metric, these results 

suggest AM communities are continuing to be affected by 

acidic deposition despite the modest long-term improvements 

observed in pH. This conclusion is consistent with results 

of acid-base chemistry modeling studies in the park (Cosby 

and others, 2006; Sullivan and others, 2008) that predict that 

recovery of park streams to acidic deposition will lag behind 

improvements in air quality, in part, because of concomitant 

declines in the base cations (calcium, magnesium, sodium, and 

potassium) that neutralize acidity. 

Bulger and others (1995) reported a positive relation in 

the park between fish richness and minimum ANC values, 

and attributed this to the loss of fish species in siliciclastic 

watersheds. Although air quality has improved over the study 

period (Fioletov and others, 2011), species richness was lower 

in siliciclastic watersheds (fig. 31) and did not increase over 

time at these sites (fig. 32). Furthermore, abundances of YOY 

and adult (age 1+) brook trout were substantially lower in 

siliciclastic watersheds than in others (figs. 33 and 35). These 

results, therefore, are consistent with Bulger and others (1995) 

and further indicate that improvements in air quality have not 

substantially influenced species richness or brook trout populations in acid-sensitive watersheds. The observed increases 

in species richness in granitic and basaltic watersheds could 

indicate that these sites (1) are recovering from earlier losses 

of species richness, (2) support more suitable habitat for 

potentially immigrating fish species, and (or) (3) are connected 

to source populations that promote colonization dynamics 

within stream networks.
Relations Among Environmental Factors, Water Quality, Aquatic Macroinvertebrates, and Fish 61

Water Temperature

Park-wide increases in mean, median, and maximum 

water temperatures were detected for multiple time periods 

over the last 30 years (fig. 10), and the magnitude of these 

trends was substantial and probably biologically meaningful. The average annual park-wide increase in mean water 

temperature was 0.04 °C, which equates to about 1.2 °C over 

the last 30 years. Short-term trends generally coincided with 

long-term trends, but were more variable. Increases in water 

temperature generally reflected increases in air temperature 

over the same time period, suggesting that 

water temperatures are responding to warmer 

air temperatures. The cause for warming air 

temperature over the last 30 years, however, 

is not clear, although observed patterns suggest that global (climate change) or regional 

(North Atlantic oscillation, land-use change) 

factors seem more likely than local (gypsy 

moth defoliation) factors as possible explanations. Irrespective of the cause, there is evidence that warming water temperatures have 

affected biological communities in the park. 

The Leuctra:P% metric seemed to track 

park-wide increases in water temperature, 

as the Leuctra:P% showed declines across 

all four time-frames evaluated and, as with 

observed increases in water temperature, 

these declines were independent of geology (fig. 24). Observed decreases in the 

Leuctra:P% metric are more consistent with 

park-wide increases in water temperatures 

than decreases in acidity. As a group, the 

stoneflies (order Plecoptera) are the most 

stenothermic of aquatic insect orders, requiring a relatively narrow range of cool or cold 

water to complete their life cycle (Sweeney and others, 1992; Fochetti and Tierno 

De Figueroa, 2008). Moreover, the genus

Leuctra is one of the most stenothermic 

genera of stonefly taxa and has been implicated as particularly vulnerable to warming 

water temperatures (Tierno de Figueroa and 

others, 2010). 

In order to examine the relation 

between water temperature and AM communities more fully, both mean annual water 

temperature and mean abundance of Leuctra 

were plotted by study year (fig. 38), showing 

that temporal patterns in water temperature 

and Leuctra abundance are essentially mirror 

images of each other. For example, during 

the most recent 5-yr period between 2005 

and 2009, water temperatures were substantially higher than average (1–3 standard 

deviations above the long-term average) 

while Leuctra abundance was substantially lower than average (1–3 standard deviations below the long-term average). 

In contrast, between 1999 and 2004, water temperatures were 

generally cooler than average, while Leuctra abundance was 

higher than average. During the early part of the monitoring 

period (1991–99), water temperatures exhibited considerably 

more annual variation than during the later part, fluctuating 

between warmer than average and cooler than average every 

1–2 years. Even during this more variable period, however, 

the abundance of Leuctra mirrored water temperatures closely 

(fig. 38). The relation between water temperature and the 

North District (N=4 sites)

Central District (N= 5 sites)

South District (N = 5 sites)

EXPLANATION

North District (N=5 sites)

EXPLANATION

Central District (N= 6 sites)

South District (N = 2 sites)

199119921993199419951996199719981999

200020012002

2003

2004

2005

2006

20072008

2009

Year

Z-score of mean annual temperature Z-score of mean abundance of Leuctra

–3

–2

–1

0

1

2

3

–4

–3

–2

–1

0

1

2

3

Figure 38. Z-scores of mean annual temperature (°C) and abundance of the 

genus Leuctra (Log abundance +1) for sites with records spanning the entire 

monitoring period, grouped by Shenandoah National Park district. Graphs 

show means ± 1 standard deviation.
62 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

percentage of the total community composed of the genus 

Leuctra also was evaluated (data not shown). The analysis 

demonstrated that the percentage of Leuctra in the community also tracked water temperatures, but not as closely as the 

total abundance of Leuctra, suggesting that variation in total 

AM abundance was influenced by factors other than water 

temperature.

The data just presented provide strong support for the 

hypothesis that park-wide declines in Leuctra:P% metric 

resulted from the effect of increasing water temperatures on 

the thermally sensitive and numerically abundant stonefly, 

Leuctra. Although these analyses cannot indicate whether 

recent regional increases in temperature are the result of climate change, the absolute or relative abundance of the genus 

Leuctra probably is responding to increasing water temperature. Therefore, the genus Leuctra may serve as a sensitive 

indicator of climate change. 

The fact that the genus Leuctra possesses traits that make 

it a useful indicator of both acidity and water temperature 

changes in streams poses a dilemma. Leuctra is relatively 

acid tolerant and the decline in its relative abundance may be 

indicative of improving water quality (decreasing acidity); 

however, it is thermally sensitive and thus its decline may be 

indicative of increasing water temperature. The majority of 

evidence presented herein supports the latter explanation. This 

dilemma also raises the question of whether a broader perspective on species traits might serve to improve trend assessment 

in the park. For example, there may be other taxa that, like 

Leuctra, are sensitive to increasing water temperatures but, 

unlike Leuctra, are also sensitive to acidity. An examination 

of trends in these other indicator taxa would help distinguish 

between the effects of multiple stressors, although such examination is beyond the scope of this study. Additional details on 

the concept of a species-trait approach to bioassessment are 

presented later herein. 

Changing water temperatures may affect stream fish communities through several direct and indirect pathways. Direct 

effects involve physiological tolerance limits (such as upper 

lethal water temperatures), growth rates, disease incidence 

rates, and biotic interactions, such as predation and competition rates. Indirect effects primarily involve the availability 

and quality of prey sources, such as aquatic macroinvertebrates and juvenile fish for adult brook trout and American 

eel, and the secondary consequences of changing fish densities and richness. It would be expected that the observed 

increases in water temperature would increase fish richness, 

and evidence in support of this prediction was found within 

granitic and basaltic watersheds (14-year trends; fig. 32). 

Water temperature increases also were observed in siliciclastic 

watersheds, however, where fish richness did not increase. It 

is suspected that colonization by fishes from outside the park 

boundaries is facilitated by increasing water temperatures, 

and such upstream colonization has been reported elsewhere 

(Roberts and Hitt, 2010; Hitt and Roberts, 2012). For example, 

observed increases in maximum water temperature (fig. 10) 

did not reach thermal maxima for most cool-water fishes 

commonly found in Virginia (Cherry and others, 1977). Further increases in fish richness would be expected if water temperatures continue to rise, reflecting the shift from a cold- to 

cool-water stream ecosystem and the colonization of streams 

in the park by source populations outside park boundaries. 

Further research is needed to predict asymptotic species richness given alternative water-temperature scenarios.

Increasing water temperatures also are expected to affect 

brook trout distributions and abundance through influences on 

growth rates, prey quality, competition/predation, and physiological tolerance limits. There currently is no evidence that 

increasing water temperatures are affecting population growth 

rates for brook trout in the park, because the observed mean 

water temperatures are generally less than 17 oC, the preferred 

limits for brook trout, (Cherry and others, 1977). Moreover, 

marginal increases in water temperature will increase individual growth rates and may increase survival, especially 

for YOY fishes. The observed increase in water temperature 

(approximately 1.2 oC over 30 years) was estimated using 

linear modeling techniques, and future changes in water temperature may or may not exhibit linearity; therefore, continued 

monitoring of brook trout populations is critical for Shenandoah National Park to manage this resource effectively, should 

environmental factors alter the rate of water temperature 

change over time.

Forest Defoliation

Forest defoliation by insect pests, in particular the gypsy 

moth (Lymantria dispar), could have affected stream condition 

and associated biological structure and composition over the 

monitoring period. The gypsy moth was first reported in the 

park in 1984, more than 60 percent defoliation was observed 

by 1987 in the northern section of the park, and defoliation had 

spread throughout the park by 1992, including the southern 

section. For example, 85.2 percent of the Paine Run watershed 

was defoliated by 1992 (Eshleman and others, 1998). Thus, 

gypsy moth effects were prevalent throughout the monitoring 

period and, although periodic defoliation from gypsy moth 

continues, much of the park shows marked recovery. 

Effects of forest defoliation by gypsy moth might be 

expected to cascade to streams in several ways. First, nutrient 

concentrations, especially that of dissolved nitrogen, might 

be expected to increase in streams, either from direct inputs 

through gypsy moth frass, or through reduced uptake of soil 

nutrients by the affected forest. The amount of light reaching streams also would be expected to increase as forests are 

defoliated. Because nutrients and light are naturally limited 

resources in headwater streams, increases might be expected 

to fuel autochthonous (instream) primary production by benthic algae. At the same time, allochthonous inputs from leaf 

material of the surrounding forest, the primary carbon source 

of headwater streams, probably would be reduced in defoliated 

watersheds. All of these factors would be expected to force a 

shift in trophic composition of AM consumers from shredders, 
Considerations for Future Monitoring 63

which capitalize on large leaf material that fall into the stream, 

to scrapers, which feed primarily on benthic algae (Vannote 

and others, 1980). Lastly, the increase in sunlight associated 

with insect defoliation could increase water temperatures. 

Therefore, it is possible that forest defoliation might explain 

the observed increases in water temperature and associated 

AM community responses described earlier. 

Previous research conducted in selected watersheds 

in the park indicated that gypsy moth defoliation results in 

sharp increases in NO3

–

 concentrations in streams, but these 

increases occur primarily during storm events (Eshleman 

and others, 1998; Scanlon and others, 2010). As previously 

described, it was not possible to use the available NO3

–

 data to 

conduct any meaningful analysis of NO3

-

 trends because of the 

preponderance of non-detects in the dataset and limitations of 

the modeling approach. Based on these earlier intensive studies, it is reasonable to assume that storm-event NO3

–

 concentrations have increased substantially in park streams relative 

to pre-infestation periods, but one would expect some level of 

recovery during the more recent part of the monitoring period. 

Although light levels have not been systematically monitored 

in the park, it is possible that light levels have increased relative to pre-infestation periods as well.

Aquatic macroinvertebrate trophic metrics (Shredder% 

and Scraper%) showed significant temporal trends in park 

streams, although temporal patterns do not appear to be related 

to gypsy moth induced changes in stream nutrient concentrations. Specifically, although the Shredder% metric declines 

throughout the monitoring period, the Scraper% metric also 

shows long-term declines (fig. 25). The declines in both metrics were limited to sites in siliciclastic watersheds because no 

significant changes were observed in granitic or basaltic watersheds (fig. 25). Unless an argument could be made that gypsy 

moth defoliation was more prevalent in siliciclastic watersheds 

than others, the observed patterns are not consistent with 

trophic responses to nutrient or light enhancement caused by 

gypsy moth defoliation. Previous research suggests that biotic 

uptake of nitrogen by stream communities can be substantial, 

and shifts in stream trophic structure in response to nutrient fluxes caused by forest disturbances have been observed 

(Webster and others, 1992; Bernhardt and others, 2003). 

The efficiency of stream communities in removing nitrogen, 

however, can be significantly influenced by stream volume and 

flow (Arheimer and others, 1996; Marti and others, 1997). In 

particular, the biological uptake of nitrogen has been shown 

to be relatively unimportant during storm events because of 

reduced residence time and decreased contact between nutrient-laden water and stream sediment, the site of biological 

uptake in lotic systems (Arheimer and others, 1996; Marti and 

others, 1997). Because the effects of gypsy moth defoliation 

on stream NO3

–

 concentrations were mostly observed during 

storm events (Eshleman and others, 1998; Scanlon and others, 

2010), it is likely that the elevated NO3

–

 concentrations were 

not available to headwater communities in the park, although 

effects could be substantial in larger stream reaches or lakes 

and ponds farther downstream.

As discussed earlier, the park-wide decline observed in 

the Leuctra:P% metric was strongly correlated with temporal 

increases in water temperatures. Although climate change is 

probably not implicated, water temperature generally tracked 

air temperature. Nevertheless, it is possible that forest defoliation caused by the gypsy moth may have contributed to the 

observed declines in the Leuctra:P% metric by enhancing 

water-temperature increases. In this assessment of temporal 

patterns in water temperature (fig. 38A), however, no obvious 

spatial pattern to the increasing water temperature trend that 

would indicate a gypsy moth effect was found. Specifically, 

gypsy moth defoliation began in the north section of the park 

and proceeded south over time; consequently, if defoliation was an important factor in observed water-temperature 

increases, streams in the north district would be expected to be 

warmer early in the monitoring period and streams in the south 

district would be expected to be warmer later in the monitoring period. Instead, the warmest district generally seems to be 

random in location across years (fig. 38A). 

The relation between watershed area and EPT%, Shredder%, PTV, and Intol% observed for the most recent 5-year 

time period (fig. 26) did not appear to correlate with water 

quality or water temperature. Because the patterns observed 

for all four of these AM metrics indicate that smaller watersheds showed greater declines in stream condition than larger 

ones, it is hypothesized that these short-term trends were due 

to more severe or more frequent flow disturbance events, such 

as droughts and (or) floods. A rigorous assessment of the relations between AM community structure and temporal variation 

in flow would require defining flow disturbances in a biologically relevant way, which is beyond the scope of this study. 

The data needed to conduct this type of analysis are available. 

The ability to identify and account for the effects of natural disturbances, such as droughts and floods, on biological 

metrics could substantially improve the power to detect trends 

associated with anthropogenic disturbances. 

Considerations for Future Monitoring

The current monitoring design in the park is spatially and 

temporally robust, specifically, in terms of the ability it provides to detect changes in water quality and aquatic fauna and 

to understand processes related to stream acidification. The 

current design generates the requisite data to satisfy multiple 

objectives, including those for which the monitoring program 

was not originally designed (for example, water-temperature 

trend analysis). Changes to this design, however, may allow 

for more efficient accomplishment of monitoring objectives or more complete representation of water-resource and 

aquatic-fauna conditions throughout the park than is currently 

possible. Some of these potential changes apply to the entire 

monitoring program, whereas others apply to specific program 

components.
64 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

The overall monitoring program in the park could be 

improved by better integration of the three monitored components, namely, water quality, AM, and fish. Currently, the 

program is well designed to evaluate spatial and temporal 

trends in each component, but managers are limited in drawing 

inferences regarding causal linkages, because each component 

is monitored using a different set of sites. Thus, one possible improvement to the program would be to co-locate the 

data-collection efforts for multiple components at more sites. 

Co-located monitoring stations for water-quality, AM, and fish 

sampling would improve the ability to perform analyses of 

interdependencies between these monitored components of the 

aquatic system.

Five long-term streamgages are operated within the park 

as part of the SWAS program, generating valuable continuous 

streamflow data to aid in the interpretation of water-quality 

and aquatic biota data. Most watersheds in the park do not 

have a streamgage, however, and although extrapolation methods can be used to represent general hydrologic conditions, the 

most effective data analysis approaches rely upon streamflow 

data collected within the watershed of interest. Additionally, 

few data are available to describe the location and timing 

of natural disturbance effects, such as those of drought (and 

consequent periods of no flow) or floods in the park. It is 

hypothesized that a large fraction of short-term (5-year) temporal variation in AM community structure observed in these 

analyses was due to the responses of AM communities to flow 

disturbance events. The ability to account for natural disturbance events may provide more power to detect change caused 

by anthropogenic stressors, which is the primary purpose of 

the program. Therefore, additional streamgages in currently 

ungaged areas would benefit all of the aquatic monitoring 

programs within the park.

Water-Quality Monitoring Program

One change to the current water-quality monitoring program that would result in more efficient data collection, while 

still accomplishing the objective of documenting long-term 

trends in water quality, is to reduce the frequency of sample 

collection at the sites currently sampled on a weekly basis. 

As noted earlier in the temporal-trends discussion, the highest frequency of sampling events that can be used for trend 

analysis is a monthly time step. Therefore, weekly samples are 

unnecessary for the evaluation of long-term temporal trends 

and cost savings would result if the sites sampled weekly were 

instead sampled monthly. Furthermore, trend analysis of quarterly sampling data has been shown herein (fig. 6) to generate 

results comparable to those computed from monthly sampling 

data; therefore, a further reduction in the frequency of sample 

collection is achievable, while still meeting the objective of 

documenting long-term trends. These changes could be made 

to the monitoring program to save costs, if required, or to 

generate the funding needed to support additional monitoring 

stations and spatially expand the trend-monitoring network. 

Such an expansion would be beneficial, because an expanded 

monitoring network would provide a more complete picture 

of long-term trends in water-resource conditions in the park. 

It is important to recognize, however, that weekly sampling 

may be valuable to satisfy different objectives of the monitoring effort and that the potential changes described here 

focused solely on data collection to support long-term trend 

analysis. For example, studies of contaminant flux require 

the collection of samples over the entire range of hydrologic 

conditions, and one approach to generating such a dataset is 

frequent sampling, particularly during stormflow events (for 

example, Riscassi and Scanlon, 2011). Although it is difficult 

to predict future data needs, the potential need for such studies 

should be carefully considered prior to altering the current 

monitoring program.

Aquatic Macroinvertebrate Monitoring Program

 Significant spatial and temporal trends in AM metrics 

were detected, suggesting that sampling effort (number of 

samples per site, fixed area instead of fixed count sampling), 

taxonomic level of resolution, and sampling equipment (sampling device and associated net mesh size) were sufficient to 

meet the status and trends objectives of this study. Therefore, 

continued use of these methods is expected to continue to 

provide the requisite data for evaluating trends. The taxonomic identification method used for all samples, however, is 

described as “lowest practical,” but numerous insect taxa are 

identified to the species level. For quantitative analysis of status and trends, this approach resulted in extensive taxonomic 

ambiguity at the species level, where, in almost all cases, 

taxa identified to the species level were also identified to 

higher levels of taxonomic resolution within the same sample. 

Thus, quantitative analysis of trends required rolling species 

level IDs up to the genus level to account for the variability 

in taxonomic level, thereby reducing the value of identifying some individuals to the species level. Consequently, cost 

savings potentially could be obtained by identifying taxa no 

further than genus level. Species-level identification, however, 

potentially allows park personnel or other scientists to address 

other research objectives in future studies, and such objectives should be considered prior to changing the identification 

protocol. For example, the species trait approach for bioassessment described earlier has much greater potential when considering taxa at the species level because there is considerably 

less variability in traits, including pollution tolerances (Lenat 

and Resh, 2001), within species than within higher taxonomic 

groupings. In addition, national parks are important repositories of native biodiversity and species are the fundamental 

unit of diversity. Thus, even qualitative assessments of species 

lists over time can have great value in terms of biological 

conservation. 

Park managers currently rely on a set of AM metrics 

considered to be generalized measures of stream condition 

or stress. This metric approach is the most common tool 
Considerations for Future Monitoring 65

involving AM data used for bioassessment (Jacobson, 2000; 

Lücke and Johnson, 2009). Given the scope and quality of the 

AM data collected in the park, however, it may be possible 

to develop stress-specific metrics that rely on autecological 

relations between traits of specific taxa and environmental 

gradients, such as the Leuctra:P% metric used in these analyses, that are designed to be a specific indicator of a specific 

stressor. Species-traits research is a promising area of research 

in the field of bioassessment but requires large spatial and 

temporal databases to establish pertinent linkages; the park’s 

AM (and fish) datasets offer such an opportunity. Several of 

the current metrics used should be adjusted for the stonefly 

genus Leuctra. The classification of the genus as both tolerant 

to acidic conditions and sensitive to most other stressors may 

serve more to confuse rather than elucidate inferences regarding trends. 

The panel approach to site selection used in the park currently employs approximately 26 fixed sites that are sampled 

annually and another 6 sites sampled on a rotational basis 

every 3 to 9 years. The rotational sites were initially included 

so that, over time, all fish-bearing streams would be sampled 

(J.E.B. Wofford, National Park Service, written commun., 

2011). At this point, most, if not all, fish-bearing streams have 

been sampled. From a status and trends perspective, the rotational sites do not offer much value, because they do not have 

the same level of temporal representation as the fixed sites. 

Moreover, it is clear at this point that the selected fixed sites 

are representative of the pertinent landscape characteristics 

(geology and watershed area) in the park for detecting changes 

in stream condition caused by regional and global stressors. 

Consequently, eliminating rotational sites for cost savings, or 

redistributing effort to accommodate more fixed sites might be 

appropriate at this stage of the program. 

The current suite of sites sampled for AM represents 

important landscape gradients. Small first-order headwater 

streams, however, are under-represented in the Shenandoah 

National Park AM Program. For example, first-order streams 

compose nearly 70 percent of the total stream length in the 

park (conservatively determined from blue lines on 1:24,000-

scale topographic maps), yet only 30 percent of the AM fixed 

sites are first order (Snyder and others, 2013). Only 1 of the 24 

fixed sites (those used in temporal analyses) drains an area less 

than 100 ha (site 1F308; table 9). It is clear from comparing 

results from this effort with those of Snyder and others (2013) 

that there are distinct differences in AM communities and their 

relations to landscape factors in streams having a contributing 

watershed area less than 100 ha compared to those having contributing areas greater than 100 ha. These differences are probably due to the relatively high vulnerability of small streams to 

flow disturbances, especially droughts. Consequently, increasing sampling effort in smaller watersheds would permit more 

rigorous evaluation of the interaction between natural disturbance and anthropogenic stressors. To be useful, however, 

more streamgages would be needed to monitor variation in 

streamflow throughout the park. 

Fish Monitoring Program

The park’s fish dataset provides a unique opportunity 

for understanding spatial and temporal changes in stream fish 

communities and brook-trout population dynamics. The analysis provided a synthesis of major spatial and temporal patterns 

and has several implications for the monitoring design and 

future research needs of the park. 

Long-term data were necessary to observe trends in 

fish communities and brook trout populations. Without 

exception, different trends were observed for the 5-year and 

14-year analyses, because extensive interannual variation 

in brook trout abundances dominated statistical models for 

5-year timespans, but the inclusion of additional years of 

data reduced such background variation such that the effects 

of geology and watershed area could be detected. Therefore, 

placing the highest priority for additional sampling on the sites 

with the longest existing records would permit future analysis 

of long-term trends.

Different trends were observed for YOY and adult brook 

trout, and a combined analysis may have obscured these 

fundamental differences. Therefore, continuing the sampling 

protocol, whereby YOYs are enumerated and distinguished 

from adult fish in the dataset, would permit such differentiation in the future. Separation of YOY and adult brook trout is a 

critical component of the park’s fish monitoring program.

The 3-pass sampling methodology provides a more 

accurate sample than would be possible using fewer passes. 

Preliminary analyses (not shown) indicated that 1-pass data 

missed present species in 89 (15 percent) of the 578 sampling 

events. On one occasion (2F093 in 1999), 4 species were captured during the second and third passes that were not detected 

on the first pass. The stream fish community component of 

monitoring in the park is crucial for understanding and managing the aquatic resources of the park, and given expected 

increases in fish richness over time, the 3-pass sampling 

methodology is appropriate and should not be replaced with 

1-pass samples.

The current sample reach length (100 m) may affect the 

representation of fish population and community dynamics, 

as these often are influenced by stream network-scale factors 

(Fausch and others, 2002; Hitt and Angermeier, 2011). The 

number of mesohabitat sequences included in a given 100-m 

sample reach may vary by stream size, and this could affect 

the opportunity to sample low-abundance or rare taxa. For 

example, it is expected that fewer individual pools, riffles, 

and runs would be encountered within an equal distance as 

stream volume increases. Additional sampling to quantify fish 

community composition surrounding sample sites in a subset 

of the park sample sites may aid in defining any effect of the 

reach length on sampling results. 

Fish community and brook trout responses to changing 

water temperature and acidification merit continued research. 

First, climate-change models predict increasing maximum 

surface-water temperatures, and long-term park data support this prediction. The implications for fish communities 
66 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

and populations, however, are uncertain. Additional research 

would be necessary to develop spatially explicit models 

relating air and water temperatures. This work would provide 

inferences about the role of groundwater and the relative vulnerability of stream segments to increasing air temperatures. 

In turn, this would also enable spatially explicit predictions 

for brook trout occupancy in the future. Second, population growth rate is an important parameter for assessing fish 

population responses to any environmental factor (including 

acid deposition and temperature change), and such assessment 

requires sequential years of data. Given the observed spatially 

synchronous population growth rates park-wide, it is expected 

that reduction of the currently implemented biannual sampling 

to an annual sampling protocol would offer important benefits 

for trend detection without loss of information as a result of 

sampling fewer locations each year.

Summary

Three disparate datasets were compiled into a comprehensive database to support the analyses conducted by this 

study, to support future efforts to evaluate aquatic resources in 

Shenandoah National Park, and to provide a working database 

consistent with National Park Service protocols for the storage 

of data collected in the future. The resulting database, which 

was created within the NPSTORET framework, contains measurements of water quality, stream habitat, aquatic macroinvertebrates (AM), and fish, from over 600 sites and covering a 

time span of over 30 years.

The assessment of stream-water quality in the park 

presented herein found a conflicting picture of acidity trends. 

Long-term trends in acid-neutralizing capacity (ANC) and 

pH suggest some improvement (decreasing acidity), but 

short-term trends suggest increasing acidity. In addition, the 

long-term increase in pH occurred park-wide, although the 

increases were weak in siliciclastic watersheds. In contrast, 

the long-term increases in ANC were largely limited to sites 

in basaltic drainages (figs. 8 and 9). SO4

2– concentration, a 

constituent of acidic deposition, indicated long-term improvements in basaltic watersheds (decreases in SO4

2– concentrations), increasing concentrations in granitic watersheds, but 

no trend in siliciclastic or basaltic watersheds (fig. 9). Overall, 

this analysis demonstrates that watersheds with higher mean 

elevations, lower land-surface gradients, and greater proportions of basaltic and carbonate geology, are least affected by 

acidification and tend to be improving over time. The watersheds that are affected by acidification and are experiencing 

continued degradation in water quality have greater proportions of siliciclastic and granitic geology, with smaller watershed areas and higher minimum watershed elevations.

In terms of water temperature, park-wide increases 

in mean, median, and maximum water temperatures were 

detected over the last 20 years (fig. 10), and the magnitude 

of these trends was substantial and probably biologically 

meaningful. The average annual increase in mean water temperature park-wide was 0.04 °C, which equates to about 1.2 °C 

over the last 30 years. Short-term trends generally coincided 

with long-term trends, but were more variable. Although 

increases in water temperatures generally tracked increasing 

air temperatures, additional analyses of longer-term (more 

than 80-year) regional air-temperature data showed that the 

observed increases in air temperature have occurred at other 

times in the past; therefore, it is not clear what mechanism is 

responsible for the observed increases over the last 20 years.

Aquatic macroinvertebrate (AM) community metrics 

showed strong evidence of temporal trends in stream condition. Although the overall condition of streams in the park 

are considered by most measures to be relatively healthy, 

there is evidence from spatial and temporal trend assessments that streams in siliciclastic watersheds have been, and 

continue to be, affected by acidic deposition, despite regional 

improvements in air quality. Although climate change is not 

conclusively implicated, streams in the park have warmed 

significantly over the last 20 years, and there is evidence that 

AM communities have responded to the warming trend. The 

correlation provides some assurance that AM monitoring 

protocols and metrics are sufficiently sensitive to detect future 

climate-change effects. There is no clear evidence that trends 

in AM communities are responding, or have responded, to 

gypsy moth defoliation.

The quality of the AM data collected over the 20-year 

period is high (consistent sampling methods employed and 

taxonomic ID conducted by the same personnel); therefore, 

confidence in the reported patterns and relations is high. Some 

of the limitations of the collected data and associated exploratory analyses that were employed are worth noting. First, data 

for water quality, AM, and fish were derived from different 

sites, and therefore, relations could not be directly assessed. 

Thus, although relatively strong inferences were derived from 

comparing common responses of these different responses to 

spatial and temporal gradients, better integration of sampling 

could improve the strength of trend analyses. Second, the 

selection of the four time periods selected for evaluation in 

this study essentially was arbitrary, based largely on periods 

of time for which data for stream water quality, AM, and fish 

were available. No effort was made to objectively define the 

boundaries (years) of a consistent trend. For example, only 

short-term trends for the most recent 5-year period between 

2005 and 2009 were evaluated, although it was clear from 

plotting the raw metric data by time (not shown) that there 

were numerous short-term fluctuations in the abundance 

and diversity of AM communities dispersed throughout the 

20-year monitoring period. It is hypothesized that a large 

fraction of the short-term variation in AM community structure may be due to the effects of natural disturbances, such as 

floods and droughts.

This analysis provided several key inferences regarding 

fish community and brook-trout population dynamics within 

the park. Increasing trends in species richness from 1996 to 

2010 were detected, and these increases were greatest in the 
References Cited 67

largest streams. Increases in richness may partially be due to 

increased water temperatures, but future research would be 

necessary to test this hypothesis. These increases in richness 

were detected only in watersheds underlain by granitic and 

basaltic geology; no trend in richness was detected in siliciclastic watersheds, which also had the lowest richness values. 

The low richness values and lack of improvement over time in 

siliciclastic watersheds is further evidence that recovery from 

acidification is not yet occurring in these watersheds. 

 Brook trout metrics provided further evidence of the 

effects of acidification in siliciclastic watersheds, because 

mean YOY and adult (age 1+) abundances were consistently 

lower in these watersheds than in others. Long-term (15-year) 

increases were detected in adult abundances in the siliciclastic 

watersheds, although YOY abundances showed a declining 

trend during this same period. These conflicting trends in YOY 

abundance and adult (age 1+) abundances demonstrate the 

importance of partitioning YOY and adult fish data in surveys. 

Increasing water temperatures apparently did not affect brooktrout population growth rates from 1996 to 2010, but continued monitoring would be necessary to evaluate future watertemperature changes and brook trout population responses.

Finally, the monitoring design in the park was found to be 

generating a spatially and temporally robust dataset, specifically in terms of the ability it provides to detect long-term 

changes in water quality and aquatic fauna, and to understand 

processes related to stream acidification. This dataset was 

found to be sufficient to satisfy multiple objectives, including 

objectives or data analysis for which the monitoring program 

was not originally designed (for example, water-temperature 

trend analysis). Multiple changes to this design that may allow 

for more efficient accomplishment of monitoring objectives 

or more complete representation of water-resources and 

aquatic fauna conditions throughout the park were identified. For example, results of the water-quality trend analysis 

demonstrate that the quarterly sampling scheme generates the 

requisite data for long-term trend monitoring, and therefore, 

savings may be achieved by reducing the sampling frequency 

for trend monitoring from weekly or monthly to quarterly. 

Discontinuities in the sampling locations for water-quality, 

AM, and fish monitoring programs, however, presented 

challenges for holistic analysis of relations between the three 

components, and future co-location of monitoring stations 

would enable more detailed analysis of the relations between 

water-quality, AM, and fish. Lastly, the scarcity of streamflow 

data collected within the park presents challenges in evaluating relations between hydrology, water-quality, AM, and 

fish—a challenge that could be overcome with the operation of 

additional streamgages.

References Cited

Arheimer, B., Andersson, L., and Lepistö, A., 1996, Variation 

of nitrogen concentration in forest streams—Influences of 

flow, seasonality and catchment characteristics: Journal of 

Hydrology, v. 179, no. 1–4, p. 281–304.

Atkinson, J.B., 2002, Shenandoah National Park fisheries 

monitoring protocol: Luray, Va., Shenandoah National Park, 

67 p.

Barbour, M.T., Gerritsen, J., Snyder, B.D., and Stribling, J.B., 

1999, Rapid bioassessment protocols for use in streams and 

wadeable rivers—Periphyton, benthic macroinvertebrates, 

and fish: Washington, D.C., U.S. Environmental Protection 

Agency, Office of Water, 339 p.

Berg, L.Y., and Moore, R.B., 1941, Forest cover types of 

Shenandoah National Park, Virginia: U.S. Department of 

the Interior, National Park Service, Region One, 132 p.

Bernhardt, E.S., Likens, G.E., Buso, D.C., and Driscoll, 

C.T., 2003, In-stream uptake dampens effects of major 

forest disturbance on watershed nitrogen export: Proceedings of the National Academy of Sciences, v. 100, no. 18, 

p. 10304–10308.

Bowman, M.F., Somers, K.M., Reid, R.A., and Scott, L.D., 

2006, Temporal response of stream benthic macroinvertebrate communities to the synergistic effects of anthropogenic acidification and natural drought events: Freshwater 

Biology, v. 51, no. 4, p. 768–782.

Bricker, O.P., and Rice, K.C., 1989, Acidic deposition to 

streams—A geology-based method predicts their sensitivity: Environmental Science & Technology, v. 23, no. 4, 

p. 379–385.

Brown, B.L., Swan, C.M., Auerbach, D.A., Campbell Grant, 

E.H., Hitt, N.P., Maloney, K.O., and Patrick, C., 2011, 

Metacommunity theory as a multispecies, multiscale 

framework for studying the influence of river network 

structure on riverine communities and ecosystems: Journal 

of the North American Benthological Society, v. 30, no. 1, 

p. 310–327.

Bulger, A.J., Dolloff, C.A., Cosby, B.J., Eshleman, K.N., 

Webb, J.R., and Galloway, J.N., 1995, The “Shenandoah 

National Park: Fish in sensitive habitats” (SNP: FISH) project. An integrated assessment of fish community responses 

to stream acidification: Water, Air, & Soil Pollution, v. 85, 

no. 2, p. 309–314.

Burns, D.A., Lynch, J.A., Cosby, B.J., Fenn, M.E., Baron, 

J.S., U.S. EPA Clean Air Markets Division, 2011, National 

Acid Precipitation Assessment Program Report to Congress 2011—An Integrated Assessment: Washington, D.C., 

National Science and Technology Council, 114 p.
68 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Burton, J., and Gerritsen, J., 2003, A stream condition index 

for Virginia non-coastal streams: Owings Mills, Md., Tetra 

Tech, Inc., 163 p.

Cherry, D.S., Dickson, K.L., Cairns, J., Jr., and Stauffer, J.R., 

1977, Preferred, avoided, and lethal temperatures of fish 

during rising temperature conditions: Journal of the Fisheries Research Board of Canada, v. 34, no. 2, p. 239–246.

Clarke, A., Mac Nally, R., Bond, N., and Lake, P.S., 2008, 

Macroinvertebrate diversity in headwater streams—A 

review: Freshwater Biology, v. 53, no. 9, p. 1707–1721.

Clarkson, R.B., 1964, Tumult on the mountains: Lumbering in 

West Virginia, 1770–1920: Parsons, W. Va., McClain Printing Co., 410 p.

Cosby, B.J., Ryan, P.F., Webb, J.R., Hornberger, G.M., and 

Galloway, J.N., 1991, Mountains of western Virginia, chap. 

10 of Charles, D.F., and Christie, S., eds., Acidic deposition and aquatic ecosystems: New York, Springer-Verlag, 

p. 297–318.

Cosby, B.J., Webb, J.R., Galloway, J.N., and Deviney, F.A., 

2006, Acidic deposition impacts of natural resources in 

Shenandoah National Park: National Park Service Technical 

Report NPS/NER/NRTR-2006/066, 35 p.

Cuffney, T.F., Bilger, M.D., and Haigler, A.M., 2007, Ambiguous taxa—Effects on the characterization and interpretation 

of invertebrate assemblages: Journal of the North American 

Benthological Society, v. 26, no. 2, p. 286–307.

Cuffney, T.F., and Brightbill, R.A., 2011, User’s manual for 

the National Water-Quality Assessment Program Invertebrate Data Analysis System (IDAS) software, version 5: 

U.S. Geological Survey Techniques and Methods book 7, 

chap. C4, 126 p.

Deviney, F.A., Jr., Brown, D.E., and Rice, K.C., 2012, Evaluation of Bayesian estimation of a hidden continuous-time 

Markov chain model with application to threshold violation in water-quality indicators: Journal of Environmental 

Informatics, v. 19, no. 2, p. 70–78. (Also available at http://

dx.doi.org/10.3808/jei.201200210.) 

Deviney, F.A., Jr., Rice, K.C., and Hornberger, G.M., 2006, 

Time series and recurrence interval models to predict 

the vulnerability of streams to episodic acidification in 

Shenandoah National Park, Virginia: Water Resources 

Research, v. 42, W09405, (Also available at http://dx.doi.

org/10.1029/2005WR004740.)

DeVivo, M.S., 1986, The deforestation of western North Carolina—1900 to 1920: Pioneer America Society Transactions, 

v. 9, p. 89–94.

Eshleman, K.N., Morgan, R.P., II, Webb, J.R., Deviney, F.A., 

and Galloway, J.N., 1998, Temporal patterns of nitrogen 

leakage from mid-Appalachian forested watersheds: Role of 

insect defoliation: Water Resources Research, v. 34, no. 8, 

p. 2005–2016.

Fausch, K.D., Torgersen, C.E., Baxter, C.V., and Li, H.W., 

2002, Landscapes to riverscapes: Bridging the gap between 

research and conservation of stream fishes: BioScience, 

v. 52, no. 6, p. 483–498.

Fenneman, N.M., and Johnson, D.W., 1946, Physical divisions 

of the United States (Map): Washington, D.C., U.S. Geological Survey.

Fioletov, V.E., McLinden, C.A., Krotkov, N., Moran, M.D., 

and Yang, K., 2011, Estimation of SO2 emissions using 

OMI retrievals: Geophysical Research Letters, v. 38, no. 21, 

p. L21811.

Fjellheim, A., and Raddum, G.G., 1990, Acid precipitation—

Biological monitoring of streams and lakes: Science of The 

Total Environment, v. 96, no. 1–2, p. 57–66.

Fochetti, R., and Tierno de Figueroa, J.M., 2008, Global 

diversity of stoneflies (Plecoptera; Insecta) in freshwater: 

Hydrobiologia, v. 595, no. 1, p. 365–377.

Foster, D., Motzkin, G., O’Keefe, J., Boose, E., Orwig, D., 

Fuller, J., and Hall, B., 2004, The environmental and human 

history of New England, in Foster, D.R., and Aber, J.D., 

eds., Forests in Time: The environmental consequences of 

1,000 years of change in New England: New Haven, Conn., 

Yale University Press, p. 43–101.

Griffith, M.B., Perry, S.A., and Perry, W.B., 1995, Macroinvertebrate communities in headwater streams affected by 

acidic precipitation in the central Appalachians: Journal of 

Environmental Quality, v. 24, no. 2, p. 233–238.

Haggerty, S.M., Batzer, D.P., and Jackson, C.R., 2002, Macroinvertebrate assemblages in perennial headwater streams of 

the Coastal Mountain range of Washington, U.S.A.: Hydrobiologia, v. 479, no. 1, p. 143–154.

Heino, J., 2005, Functional biodiversity of macroinvertebrate assemblages along major ecological gradients of 

boreal headwater streams: Freshwater Biology, v. 50, no. 9, 

p. 1578–1587.

Helsel, D.R., and Hirsch, R.M., 2002, Statistical methods in 

water resources: U.S. Geological Survey Techniques of 

Water-Resources Investigations, book 4, chap. A3, 523 p.

Hem, J.D., 1985, Study and interpretation of the chemical 

characteristics of natural water (3d ed.): U.S. Geological 

Survey Water-Supply Paper 2254, 263 p.
References Cited 69

Herlihy, A.T., Kaufmann, P.R., Church, M.R., Wigington, P.J., 

Jr., Webb, J.R., and Sale, M.J., 1993, The effects of acidic 

deposition on streams in the Appalachian Mountain and 

Piedmont Region of the Mid-Atlantic United States: Water 

Resources Research, v. 29, no. 8, p. 2687–2703.

Hirsch, R.M., and Slack, J.R., 1984, A nonparametric trend 

test for seasonal data with serial dependence: Water 

Resources Research, v. 20, no. 1, p. 727–732.

Hirsch, R.M., Slack, J.R., and Smith, R.A., 1982, Techniques 

of trend analysis for monthly water quality data: Water 

Resources Research, v. 18, no. 1, p. 107–121, (Also available at http://dx.doi.org/10.1029/WR018i001p00107.)

Hitt, N.P., and Angermeier, P.L., 2011, Fish community 

and bioassessment responses to stream network position: 

Journal of the North American Benthological Society, v. 30, 

no. 1, p. 296–309.

Hitt, N.P., Eyler, S., and Wofford, J.E.B., 2012, Dam removal 

increases American Eel abundance in distant headwater 

streams: Transactions of the American Fisheries Society, 

v. 141, no. 5, p. 1171–1179.

Hitt, N.P., and Roberts, J.H., 2012, Hierarchical spatial 

structure of stream fish colonization and extinction: Oikos, 

v. 121, no. 1, p. 127–137.

Homer, Collin, Huang, Chengquan, Yang, Limin, Wylie, Bruce 

and Coan, Michael, 2004, Development of a 2001 National 

Landcover Database for the United States: Photogrammetric 

Engineering and Remote Sensing, v. 70, no. 7, p. 829–840.

Jacobson, P.T., 2000, Evaluation of multi-metric bioassessment as an approach for assessing impacts of entrainment 

and impingement under Section 316(b) of the Clean Water 

Act: Environmental Science & Policy, v. 3, supplement 1, 

p. 107–115.

Jenkins, R.E., and Burkhead N.M., 1994, The freshwater 

fishes of Virginia: Bethesda, M.D., American Fisheries 

Society, 1079 p.

Kalnay, E., and Cai, M., 2003, Impact of urbanization and 

land-use change on climate: Nature, v. 423, p. 528–531.

Kaushal, S.S., Likens, G.E., Jaworski, N.A., Pace, M.L., 

Sides, A.M., Seekell, D., Belt, K.T., Secor, D.H., and Wingate, R.L., 2010, Rising stream and river temperatures in the 

United States: Frontiers in Ecology and the Environment, 

v. 8, no. 9, p. 461–466.

Kendall, M.G., 1975, Rank correlation methods: Oxford, U.K., 

Griffin, 202 p. 

Kimmel, W.G., Murphey, D.J., Sharpe, W.E., and DeWalle, 

D.R., 1985, Macroinvertebrate community structure and 

detritus processing rates in two southwestern Pennsylvania 

streams acidified by atmospheric deposition: Hydrobiologia, 

v. 124, no. 2, p. 97–102.

Klemm, D.J., Blocksom, K.A., Thoeny, W.T., Fulk, F.A., 

Herlihy, A.T., Kaufmann, P.R., and Cormier, S.M., 2002, 

Methods development and use of macroinvertebrates as 

indicators of ecological conditions for streams in the MidAtlantic Highlands Region: Environmental monitoring and 

assessment, v. 78, no. 2, p. 169–212. 

Lenat, D.R., 1993, A biotic index for the southeastern United 

States: Derivation and list of tolerance values, with criteria 

for assigning water-quality ratings: Journal of the North 

American Benthological Society, v. 12, no. 3, p. 279–290.

Lenat, D.R., and Resh, V.H., 2001, Taxonomy and stream 

ecology—The benefits of genus- and species-level identifications: Journal of the North American Benthological 

Society, v. 20, no. 2, p. 287–298.

Likens, G.E., 2010, River ecosystem ecology—A global perspective: San Diego, Calif., Academic Press, 425 p.

Lücke, J.D., and Johnson, R.K., 2009, Detection of ecological change in stream macroinvertebrate assemblages using 

single metric, multimetric or multivariate approaches: 

Ecological Indicators, v. 9, no. 4, p. 659–669.

Mann, H.B., 1945, Nonparametric tests against trend: Econometrica, v. 13, no. 3, p. 245–259.

Martí, E., Grimm, N.B., and Fisher, S.G., 1997, Pre- and 

post-flood retention efficiency of nitrogen in a Sonoran 

desert stream: Journal of the North American Benthological 

Society, v. 16, no. 4, p. 805–819.

McCune, B., Grace, J.B., and Urban, D.L., 2002, Analysis 

of ecological communities: Gleneden Beach, Oreg., MjM 

Software Design, 304 p.

Merritt, R.W., and Cummins, K.W., 1996, An introduction 

to the aquatic insects of North America: Dubuque, Iowa, 

Kendall Hunt, 862 p.

Morgan, B.A., Eaton, L.S., and Wieczorek, G.F., 2003, Pleistocene and Holocene colluvial fans and terraces in the Blue 

Ridge region of Shenandoah National Park, Virginia: U.S. 

Geological Survey Open-File Report 03-410, 25 p.

National Oceanic and Atmospheric Administration, 2012, 

Climate Data Online Database: National Climatic Data 

Center, accessed March 12, 2012, at http://www.ncdc.noaa.

gov/cdo-web/.

National Park Service, 1998, Resource Management Plan: 

Luray, Va., Shenandoah National Park. 
70 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

National Park Service, 2008, Natural resource fact sheet: 

Luray, Va., Aquatic macroinvertebrate monitoring program. 

National Park Service, 2012, Weather—Shenandoah National 

Park: Accessed February 15, 2012, at http://www.nps.gov/

shen/naturescience/weather.htm.

National Park Service, 2013, Fish—Shenandoah National 

Park: Accessed April 26, 2013, at http://www.nps.gov/shen/

naturescience/fish.htm

Rice, K.C., Deviney, F.A., and Olson, G., 2007, Acid rain 

in Shenandoah National Park, Virginia: U.S. Geological 

Survey Fact Sheet 2007–3057, 4 p. (Also available at http://

pubs.usgs.gov/fs/2007/3057/).

Riscassi, A.L., and Scanlon, T.M., 2011, Controls on stream 

water dissolved mercury in three mid-Appalachian forested 

headwater catchments: Water Resources Research, v. 47, 

no. 12.

Roberts, J.H., and Hitt, N.P., 2010, Longitudinal structure in 

temperate stream fish communities—Evaluating conceptual 

models with temporal data: American Fisheries Society 

Symposium, v. 73, p. 281–299.

Roghair, C.N., and Dolloff, C.A., 2005, Brook Trout movement during and after recolonization of a naturally defaunated stream reach: North American Journal of Fisheries 

Management, v. 25, no. 3, p. 777–784.

Rosenberg, D.M., and Resh, V.H., 1993, Freshwater biomonitoring and benthic macroinvertebrates: New York, Chapman 

& Hall, 516 p.

Scanlon, T.M., Ingram, S.M., and Riscassi, A.L., 2010, Terrestrial and in-stream influences on the spatial variability 

of nitrate in a forested headwater catchment: Journal of 

Geophysical Research, v. 115, no. G2, p. G02022.

Schertz, T.L., Alexander, R.B, and Ohe, D.J., 1991, The computer program EStimate TREND (ESTREND), a system for 

the detection of trends in water-quality data: U.S. Geological Survey Water-Resources Investigations Report 91–4040, 

63 p.

Shaw, P.J.A., 2003, Multivariate statistics for the environmental sciences: New York, Oxford University Press, 233 p.

Shreve, R.L., 1966, Statistical law of stream numbers: The 

Journal of Geology, v. 74, no. 1, p. 17–37.

Slack, J.R., Lorenz, D.L., and others, 2003, USGS library for 

S-PLUS for Windows—Release 2.1: U.S. Geological Survey Open-File Report 03–357

Snyder, C.D., and Johnson, Z.B., 2006, Macroinvertebrate 

assemblage recovery following a catastrophic flood and 

debris flows in an Appalachian mountain stream: Journal 

of the North American Benthological Society, v. 25, no. 4, 

p. 825–840.

Snyder, C.D., Webb, J.R., Young, J.A., and Johnson, Z.B., 

2013, Significance of headwater streams and perennial 

springs in ecological monitoring in Shenandoah National 

Park: U.S. Geological Survey Open-File Report 2013–1178, 

46 p.

Snyder, C.D., Young, J.A., and Stout, B.M., III, 2006, Aquatic 

habitats of Canaan Valley, West Virginia: Diversity and 

environmental threats: Northeastern Naturalist, v. 13, no. 3, 

p. 333–352.

Sokal, R.R., and Rohlf, F.J., 1995, Biometry (3d ed.): New 

York, W.H. Freeman and Company, 88 p.

Soucek, D.J., 2007, Sodium sulfate impacts feeding, specific 

dynamic action, and growth rate in the freshwater bivalve 

Corbicula fluminea: Aquatic toxicology, v. 83, no. 4, 

p. 315–322.

Southworth, Scott, Aleinikoff, J.N., Bailey, C.M., Burton, 

W.C., Crider, E.A., Hackley, P.C., Smoot, J.P., and Tollo, 

R.P., 2009, Geologic map of the Shenandoah National 

Park region Virginia: U.S. Geological Survey Open-File 

Report 2009–1153, 96 p., 1 plate, scale 1:100,000.

Strahler, A.N., 1952, Hypsometric (area-altitude) analysis of 

erosional topography: Geological Society of America Bulletin, v. 63, no. 11, p. 1117–1142.

Sullivan, T.J., Cosby, B.J., Webb, J.R., Dennis, R.L., Bulger, 

A.J., and Deviney, F.A., Jr., 2008, Streamwater acid-base 

chemistry and critical loads of atmospheric sulfur deposition in Shenandoah National Park, Virginia: Environmental 

Monitoring and Assessment, v. 137, no. 1–3, p. 85–99.

Sweeney, B.W., Jackson, J.K., Newbold, J.D., and Funk, D.H., 

1992, Climate change and the life histories and biogeography of aquatic insects in eastern North America, in Firth, P., 

and Fisher, S.G., eds., Global climate change and freshwater 

ecosystems: New York, Springer-Verlag Press, 321 p.

Tierno de Figueroa, J., López-Rodríguez, M., Lorenz, A., 

Graf, W., Schmidt-Kloiber, A., and Hering, D., 2010, Vulnerable taxa of European Plecoptera (Insecta) in the context 

of climate change: Biodiversity and Conservation, v. 19, 

no. 5, p. 1269–1277.

U.S. Department of the Interior, 1932, Final report of the 

Southern Appalachian National Park Commission to the 

Secretary of the Interior: Washington, D.C., United States 

Government Printing Office, 39 p.
References Cited 71

Vana-Miller, D.L., and Weeks, D.P., 2004, Shenandoah 

National Park, Virginia, water resources scoping report: 

Technical Report NPS/NRWRS/NRTR-2004/320, 138 p.

Vannote, R.L., Minshall, G.W., Cummins, K.W., Sedell, J.R., 

and Cushing, C.E., 1980, The river continuum concept: 

Canadian Journal of Fisheries and Aquatic Sciences, v. 37, 

no. 1, p. 130–137.

Venables, W.N., and Ripley, B.D., 2000, S Programming (1st 

ed.): New York, Springer-Verlag, 266 p.

Virginia Division of Mineral Resources, 2003, Digital representation of the 1993 geologic map of Virginia: Virginia 

Division of Mineral Resources Publication 174 [CD-ROM; 

2003, December 31]. [Adapted from Virginia Division of 

Mineral Resources, 1993, Geologic map of Virginia: Virginia Division of Mineral Resources, scale 1:500,000.]

Visbeck, M.H., Hurrell, J.W., Polvani, L., and Cullen, H.M., 

2001, The North Atlantic Oscillation: Past, present, and 

future: Proceedings of the National Academy of Sciences, 

v. 98, no. 23, p. 12876–12877.

Webster, J.R., Golladay, S.W., Benfield, E.F., Myer, J.L., 

Swank, W.T., and Wallace, J.B., 1992, Catchment disturbance and stream response—An overview of stream 

research at Coweeta Hydrologic Laboratory, chap. 15, in

Boon, P.J., Calow, P., and Petts, G.E., eds., River conservation and management: Chichester, U.K., John Wiley and 

Sons, Ltd., 470 p.

Wohl, D.L., Wallace, J.B., and Meyer, J.L., 1995, Benthic 

macroinvertebrate community structure, function and 

production with respect to habitat type, reach and drainage 

basin in the southern Appalachians (U.S.A.): Freshwater 

Biology, v. 34, no. 3, p. 447–464.

Zar, J.H., 1999, Biostatistical analysis (4th ed.): Upper Saddle 

River, N.J., Prentice Hall, 663 p.

Zipper, C.E., Holtzman, G.I., Darken, P.F., Gildea, J.J., and 

Stewart, R.E., 2002, Virginia USA water quality, 1978 to 

1995: Regional interpretation: Journal of the American 

Water Resources Association, v. 38, no. 3, p. 789–802.

Appendix 1 73

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring (LTEM) Program sites between 1990 and 2009, 

and 54 sites in 1995 (data used in the retrospective analyses contained in this report).—Continued 

[Table contains taxonomic descriptions as well as taxon-specific ecological characteristics required to compute summary statistics shown in table 10. -, taxa not 

defined at this level; USGS, U.S. Geological Survey. Functional groups are as follows: CF, Collector filterer; CG, Collector-gatherer; GN, Generalist feeder; 

MP, Macrophyte piercer; PR, Predator; SC, Scraper; and SH, shredder. Habit codes are as follows: BU, Burrowers; CG, Clingers; CL, Climbers; DV, Divers; 

GN, Generalist; SK, Skaters; and SP, Sprawlers]

Taxa 

code1 Phylum Class Order Family Genus Functionalgroup2

Tolerance 

values3 Habit4

5060 Platyhelminthes Turbellaria Tricladida Planariidae - CG 8 SP

6000 Nematoda - - - - CG 8 BU

C Nemerinea - - - - PR 8 SP

D2 Annelida Oligochaeta - - - CG 8 BU

D3 Annelida Hirudinea - - - PR 8 SP

E081 Arthropoda Crustacea Isopoda Asellidae - CG 8 SP

E090 Arthropoda Crustacea Amphipoda Talitridae - CG 8 CR

E092 Arthropoda Crustacea Amphipoda Gamaridae - CG 6 CR

E0A0 Arthropoda Crustacea Decapoda Cambaridae - GN 5 GN

E10001 Arthropoda Insecta Plecoptera Pteronarcyidae Pteronarcys SH 1 CR

E10100 Arthropoda Insecta Plecoptera Peltoperlidae Peltoperla SH 1 CR

E10101 Arthropoda Insecta Plecoptera Peltoperlidae Tallaperla SH 1 CR

E10201 Arthropoda Insecta Plecoptera Nemouridae Prostoia SH 5 CR

E10210 Arthropoda Insecta Plecoptera Nemouridae Amphinemura SH 3 CR

E10301 Arthropoda Insecta Plecoptera Perlidae Paragnetina PR 2 CR

E10302 Arthropoda Insecta Plecoptera Perlidae Agnetina PR 0 CR

E10311 Arthropoda Insecta Plecoptera Perlidae Hansonoperla PR 0 CR

E10312 Arthropoda Insecta Plecoptera Perlidae Attaneuria PR 1 CR

E10313 Arthropoda Insecta Plecoptera Perlidae Acroneuria PR 2 CR

E10314 Arthropoda Insecta Plecoptera Perlidae Eccoptura PR 3 CR

E10315 Arthropoda Insecta Plecoptera Perlidae Perlesta PR 5 CR

E10401 Arthropoda Insecta Plecoptera Perlodidae Yugus PR 0 CR

E10404 Arthropoda Insecta Plecoptera Perlodidae Remenus PR 0 CR

E10406 Arthropoda Insecta Plecoptera Perlodidae Malirekus PR 1 CR

E10411 Arthropoda Insecta Plecoptera Perlodidae Isoperla PR 2 CR

E10510 Arthropoda Insecta Plecoptera Chloroperlidae Alloperla PR 1 CR

E10511 Arthropoda Insecta Plecoptera Chloroperlidae Haploperla PR 1 CR

E10512 Arthropoda Insecta Plecoptera Chloroperlidae Sweltsa PR 0 CR

E10514 Arthropoda Insecta Plecoptera Chloroperlidae Suwallia PR 1 CR

E10600 Arthropoda Insecta Plecoptera Taeniopterygidae Taeniopteryx SH 5 CR

E10700 Arthropoda Insecta Plecoptera Leuctridae Leuctra SH 0 CR

E10801 Arthropoda Insecta Plecoptera Capniidae Paracapnia SH 1 CR

E11000 Arthropoda Insecta Ephemeroptera Ephemeridae Ephemera CG 4 BU

E11301 Arthropoda Insecta Ephemeroptera Ephemerellidae Serratella CG 2 CR

E11302 Arthropoda Insecta Ephemeroptera Ephemerellidae Timpanoga CG 4 CR

Appendix 1
74 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring (LTEM) Program sites between 1990 and 2009, 

and 54 sites in 1995 (data used in the retrospective analyses contained in this report).—Continued 

[Table contains taxonomic descriptions as well as taxon-specific ecological characteristics required to compute summary statistics shown in table 10. -, taxa not 

defined at this level; USGS, U.S. Geological Survey. Functional groups are as follows: CF, Collector filterer; CG, Collector-gatherer; GN, Generalist feeder; 

MP, Macrophyte piercer; PR, Predator; SC, Scraper; and SH, shredder. Habit codes are as follows: BU, Burrowers; CG, Clingers; CL, Climbers; DV, Divers; 

GN, Generalist; SK, Skaters; and SP, Sprawlers]

Taxa 

code1 Phylum Class Order Family Genus Functionalgroup2

Tolerance 

values3 Habit4

E11303 Arthropoda Insecta Ephemeroptera Ephemerellidae Drunella SC 0 CR

E11304 Arthropoda Insecta Ephemeroptera Ephemerellidae Ephemerella CG 2 CR

E11305 Arthropoda Insecta Ephemeroptera Ephemerellidae Eurylophella CG 4 CR

E11500 Arthropoda Insecta Ephemeroptera Ameletidae Ameletus SC 1 CG

E11600 Arthropoda Insecta Ephemeroptera Leptophlebiidae Leptophlebia CG 6 SP

E11601 Arthropoda Insecta Ephemeroptera Leptophlebiidae Paraleptophlebia CG 1 CR

E11603 Arthropoda Insecta Ephemeroptera Leptophlebiidae Habrophlebia CG 1 CR

E11604 Arthropoda Insecta Ephemeroptera Leptophlebiidae Habrophlebiodes SC 1 CR

E11701 Arthropoda Insecta Ephemeroptera Baetidae Baetis CG 5 CG

E11702 Arthropoda Insecta Ephemeroptera Baetidae Callibaetis CG 10 CL

E11703 Arthropoda Insecta Ephemeroptera Baetidae Centroptilum CG 7 CG

E11900 Arthropoda Insecta Ephemeroptera Heptageniidae Stenonema SC 4 CG

E11901 Arthropoda Insecta Ephemeroptera Heptageniidae Stenacron CG 4 CG

E11902 Arthropoda Insecta Ephemeroptera Heptageniidae Epeorus CG 1 CG

E11903 Arthropoda Insecta Ephemeroptera Heptageniidae Cinygmula SC 1 CG

E11904 Arthropoda Insecta Ephemeroptera Heptageniidae Leucrocuta SC 2 CG

E11905 Arthropoda Insecta Ephemeroptera Heptageniidae Heptagenia SC 3 CG

E11907 Arthropoda Insecta Ephemeroptera Heptageniidae Rhithrogena CG 1 CG

E11B00 Arthropoda Insecta Ephemeroptera Potomanthidae Anthopotomus CG 2 BU

E11D04 Arthropoda Insecta Ephemeroptera Isonychiidae Isonychia CF 3 CG

E12000 Arthropoda Insecta Odonata Cordulegastridae Cordulegaster PR 3 BU

E12020 Arthropoda Insecta Odonata Gomphidae Progomphus PR 8 BU

E12022 Arthropoda Insecta Odonata Gomphidae Stylogomphus PR 4 BU

E12023 Arthropoda Insecta Odonata Gomphidae Gomphus PR 6 BU

E12029 Arthropoda Insecta Odonata Gomphidae Lanthus PR 1 BU

E12041 Arthropoda Insecta Odonata Aeschnidae Boyeria PR 3 CR

E12053 Arthropoda Insecta Odonata Coenagrionidae - PR 9 CR

E131 Arthropoda Insecta Hemiptera Nepidae - PR 7 CL

E13200 Arthropoda Insecta Hemiptera Gerridae Gerris PR 8 SK

E13400 Arthropoda Insecta Hemiptera Veliidae Microvelia PR 8 SK

E13A Arthropoda Insecta Hemiptera Corixidae - MP 9 DV

E14000 Arthropoda Insecta Megaloptera Sialidae Sialis PR 7 SP

E14100 Arthropoda Insecta Megaloptera Corydalidae Corydalus PR 5 CR

E14111 Arthropoda Insecta Megaloptera Corydalidae Nigronia PR 4 CR

E16012 Arthropoda Insecta Trichoptera Hydroptilidae Hydroptila MP 6 CR
Appendix 1 75

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring (LTEM) Program sites between 1990 and 2009, 

and 54 sites in 1995 (data used in the retrospective analyses contained in this report).—Continued 

[Table contains taxonomic descriptions as well as taxon-specific ecological characteristics required to compute summary statistics shown in table 10. -, taxa not 

defined at this level; USGS, U.S. Geological Survey. Functional groups are as follows: CF, Collector filterer; CG, Collector-gatherer; GN, Generalist feeder; 

MP, Macrophyte piercer; PR, Predator; SC, Scraper; and SH, shredder. Habit codes are as follows: BU, Burrowers; CG, Clingers; CL, Climbers; DV, Divers; 

GN, Generalist; SK, Skaters; and SP, Sprawlers]

Taxa 

code1 Phylum Class Order Family Genus Functionalgroup2

Tolerance 

values3 Habit4

E16013 Arthropoda Insecta Trichoptera Hydroptilidae Ochrotrichia CG 4 CG

E16015 Arthropoda Insecta Trichoptera Hydroptilidae Stactobiella SH 1 CG

E16017 Arthropoda Insecta Trichoptera Hydroptilidae Palaeagapetus SH 5 SP

E16210 Arthropoda Insecta Trichoptera Hydropsychidae Hydropsyche CF 6 CG

E16211 Arthropoda Insecta Trichoptera Hydropsychidae Cheumatopsyche CF 6 CG

E16220 Arthropoda Insecta Trichoptera Hydropsychidae Diplectrona CF 2 CG

E16221 Arthropoda Insecta Trichoptera Hydropsychidae Homoplectra CF 2 CG

E16230 Arthropoda Insecta Trichoptera Hydropsychidae Parapsyche CF 2 CG

E16300 Arthropoda Insecta Trichoptera Rhyacophilidae Rhyacophila PR 1 CR

E16400 Arthropoda Insecta Trichoptera Philopotamidae Chimarra CF 3 CG

E16410 Arthropoda Insecta Trichoptera Philopotamidae Wormaldia CF 1 CG

E16411 Arthropoda Insecta Trichoptera Philopotamidae Dolophilodes CF 1 CG

E16500 Arthropoda Insecta Trichoptera Psychomyiidae Psychomyia CG 3 CG

E16501 Arthropoda Insecta Trichoptera Psychomyiidae Lype SC 4 CG

E167 Arthropoda Insecta Trichoptera Leptoceridae - CG 4 CR

E16A00 Arthropoda Insecta Trichoptera Odontoceridae Psilotreta SC 0 CG

E16C00 Arthropoda Insecta Trichoptera Brachycentridae Micrasema SH 1 CG

E16C01 Arthropoda Insecta Trichoptera Brachycentridae Brachycentrus CF 1 CG

E16C02 Arthropoda Insecta Trichoptera Brachycentridae Adicrophlips SH 1 CG

E16E00 Arthropoda Insecta Trichoptera Lepidostomatidae Lepidostoma SH 1 CR

E16F00 Arthropoda Insecta Trichoptera Glossosomatidae Glossosoma SC 0 CG

E16F10 Arthropoda Insecta Trichoptera Glossosomatidae Agapetus SC 0 CG

E16G10 Arthropoda Insecta Trichoptera Limnephilidae Pycnopsyche SH 4 CR

E16H00 Arthropoda Insecta Trichoptera Polycentropodidae Cyrnellus CF 7 SP

E16H01 Arthropoda Insecta Trichoptera Polycentropodidae Neureclipsis CF 4 CG

E16H02 Arthropoda Insecta Trichoptera Polycentropodidae Nyctiophylax PR 3 SP

E16H03 Arthropoda Insecta Trichoptera Polycentropodidae Polycentropus PR 4 CG

E16I00 Arthropoda Insecta Trichoptera Molannidae Molanna SC 6 SP

E16J00 Arthropoda Insecta Trichoptera Sericostomatidae Fattigia SH 1 SP

E16L00 Arthropoda Insecta Trichoptera Goeridae Goera SC 0 CR

E16M00 Arthropoda Insecta Trichoptera Uenoidae Neophylax SC 2 CG

E17 Arthropoda Insecta Lepidoptera - - 8

E182 Arthropoda Insecta Coleoptera Dytiscidae - PR 6 GN

E184 Arthropoda Insecta Coleoptera Hydrophilidae - PR 7 GN

E18700 Arthropoda Insecta Coleoptera Psephenidae Psephenus SC 4 CG
76 Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring (LTEM) Program sites between 1990 and 2009, 

and 54 sites in 1995 (data used in the retrospective analyses contained in this report).—Continued 

[Table contains taxonomic descriptions as well as taxon-specific ecological characteristics required to compute summary statistics shown in table 10. -, taxa not 

defined at this level; USGS, U.S. Geological Survey. Functional groups are as follows: CF, Collector filterer; CG, Collector-gatherer; GN, Generalist feeder; 

MP, Macrophyte piercer; PR, Predator; SC, Scraper; and SH, shredder. Habit codes are as follows: BU, Burrowers; CG, Clingers; CL, Climbers; DV, Divers; 

GN, Generalist; SK, Skaters; and SP, Sprawlers]

Taxa 

code1 Phylum Class Order Family Genus Functionalgroup2

Tolerance 

values3 Habit4

E18701 Arthropoda Insecta Coleoptera Psephenidae Ectopria SC 4 CG

E18800 Arthropoda Insecta Coleoptera Dryopidae Helichus SC 5 CG

E18900 Arthropoda Insecta Coleoptera Elmidae Stenelmis SC 5 CG

E18901 Arthropoda Insecta Coleoptera Elmidae Macronychus SC 5 CG

E18902 Arthropoda Insecta Coleoptera Elmidae Ancyronyx SC 6 CG

E18904 Arthropoda Insecta Coleoptera Elmidae Optioservus SC 5 CG

E18905 Arthropoda Insecta Coleoptera Elmidae Dubiraphia SC 5 CG

E18906 Arthropoda Insecta Coleoptera Elmidae Promoresia SC 2 CG

E18907 Arthropoda Insecta Coleoptera Elmidae Oulimnius SC 2 CG

E18908 Arthropoda Insecta Coleoptera Elmidae Gonielmas SC 5 CG

E18C00 Arthropoda Insecta Coleoptera Ptilodactylidae Anchytarsus SH 4 CG

E18F0 Arthropoda Insecta Coleoptera Helophoridae Helophorus SH 8 CL

E19000 Arthropoda Insecta Diptera Blephariceridae Blepharicera SC 0 CG

E19200 Arthropoda Insecta Diptera Tipulidae Tipula SH 5 BU

E19210 Arthropoda Insecta Diptera Tipulidae Antocha CG 3 CG

E19220 Arthropoda Insecta Diptera Tipulidae Dicranota PR 3 CR

E19222 Arthropoda Insecta Diptera Tipulidae Molophilus PR 4 BU

E19230 Arthropoda Insecta Diptera Tipulidae Hexatoma PR 3 CR

E19231 Arthropoda Insecta Diptera Tipulidae Limnophila PR 4 BU

E19234 Arthropoda Insecta Diptera Tipulidae Pilaria PR 3 BU

E19242 Arthropoda Insecta Diptera Tipulidae Ormosia CG 4 BU

E19243 Arthropoda Insecta Diptera Tipulidae Erioptera CG 3 BU

E19311 Arthropoda Insecta Diptera Psychodidae Pericoma CG 10 BU

E19500 Arthropoda Insecta Diptera Dixidae Dixa CG 3 CR

E19700 Arthropoda Insecta Diptera Simuliidae Prosimulium CF 4 CG

E19720 Arthropoda Insecta Diptera Simuliidae Simulium CF 6 CG

E198 Arthropoda Insecta Diptera Chironomidae - CG 6 BU

E199 Arthropoda Insecta Diptera Ceratopogonidae - SC 6 SP

E19B10 Arthropoda Insecta Diptera Tabanidae Tabanus PR 9 BU

E19C00 Arthropoda Insecta Diptera Athericidae Atherix PR 2 CR

E19J00 Arthropoda Insecta Diptera Empididae Hemerodromia PR 6 CR

E19J01 Arthropoda Insecta Diptera Empididae Chelifera PR 6 SP

E19J02 Arthropoda Insecta Diptera Empididae Wiedemannia PR 6 CG

E19J10 Arthropoda Insecta Diptera Empididae Clinocera PR 6 CG

E19J30 Arthropoda Insecta Diptera Empididae Oreogeton PR 6 SP
Appendix 1 77

Appendix 1. Macroinvertebrate taxa collected at 24 Long Term Ecological Monitoring (LTEM) Program sites between 1990 and 2009, 

and 54 sites in 1995 (data used in the retrospective analyses contained in this report).—Continued 

[Table contains taxonomic descriptions as well as taxon-specific ecological characteristics required to compute summary statistics shown in table 10. -, taxa not 

defined at this level; USGS, U.S. Geological Survey. Functional groups are as follows: CF, Collector filterer; CG, Collector-gatherer; GN, Generalist feeder; 

MP, Macrophyte piercer; PR, Predator; SC, Scraper; and SH, shredder. Habit codes are as follows: BU, Burrowers; CG, Clingers; CL, Climbers; DV, Divers; 

GN, Generalist; SK, Skaters; and SP, Sprawlers]

Taxa 

code1 Phylum Class Order Family Genus Functionalgroup2

Tolerance 

values3 Habit4

E1A Arthropoda Insecta Collembola - - CG 7 SK

E20 Arthropoda Arachnida Acari - - PR 6 CR

F002 Mollusca Gastropoda Limnophila Planorbidae - CG 7 SP

F014 Mollusca Gastropoda Mesogastropoda Pleuroceridae - SC 2 CG

F110 Mollusca Bivalvia Veneroida Sphaeriidae - CF 8 BU

F11100 Mollusca Bivalvia Veneroida Corbiculidae Corbicula CF 6 BU

1

Taxa codes were taken from the Shenandoah National Park long-term monitoring database.

2

Functional groups were taken from the Shenandoah National Park long-term monitoring database. These characteristics were used to compute community 

metrics Shredder% and Scraper% (see table 10).

3

Tolerance values were taken from the Shenandoah National Park long-term monitoring database. These values were used to compute community pollution 

tolerance values (see table 10). 

4

Taxon-specific “habits” were taken from the Shenandoah National Park long-term monitoring database. These characteristics were used to compute the community metric Hapto% (Haptobenthos), which is the sum of crawlers and clingers (see table 10).

Prepared by:

USGS Science Publishing Network

Raleigh Publishing Service Center

3916 Sunset Ridge Road

Raleigh, NC 27607

For additional information regarding this publication, contact:

Director

USGS Virginia Water Science Center 

1730 East Parham Road

Richmond, Virginia 23228

(804) 261-2600 

email: dc_va@usgs.gov

Or visit the Virginia Water Science Center Web site at: 

http://va.water.usgs.gov/


Jastram and others—Surface-Water Quality and Aquatic Biota Data Collected in Shenandoah National Park, Virginia—Scientific Investigations Report 2013–5157