2.1.1. Data Compilation

As of 2024, the NTN has 123 sites near USGS streamgages to monitor status (yields) and trends (change in flow-normalized loads) in nitrogen, phosphorus, and suspended sediment loads in the Chesapeake Bay watershed. Load is the product of constituent concentration and streamflow that passes a point of reference at a measured unit of time, and yield is the load normalized by watershed area. Discrete water-quality data and continuous daily streamflow data from the NTN were used to analyze the status and trends of nitrogen, phosphorus, and suspended sediment in rivers and streams of the Chesapeake Bay watershed. All data compilation and analyses were completed in R (ver. 4.1.3; R Core Team, 2022a). Streamflow time series were retrieved from the National Water Information System (NWIS) using the R package “dataRetrieval” (ver. 2.7.7; Hirsch and De Cicco, 2015; U.S. Geological Survey, 2023). Continuous (15-minute internal logging) data were used to calculate a daily mean of streamflow for each USGS streamgage. The Chesapeake Bay Program nontidal water quality monitoring program upholds high standards for NTN data quality; thus, a water year (WY) was discarded if it contained days with missing streamflow measurements that could not be remedied. Discrete nutrients (total nitrogen and total phosphorous) and total suspended sediment data were obtained from the Water Quality Portal (WQP; U.S. Geological Survey and U.S. Environmental Protection Agency, 2021). After rigorous quality-assurance checks to remove outliers and duplicates, discrete water-quality data from each site (typically 20 samples per year) were paired to the corresponding daily mean streamflows to estimate daily concentration before flow-normalization was applied for trend analysis.

2.1.2. Analysis

For nutrients and sediment, status was defined as the 10-year (WY 2011–20) average yield delivered to the Chesapeake Bay. Trends were defined as the change in flow-normalized load for a given trend interval. Trends for nitrogen, phosphorus, and suspended sediment at eligible NTN sites were estimated for two trend intervals: 1985–2020 (long term) and 2011–20 (short term). For status and the short-term trend interval, 89 total nitrogen sites, 70 total phosphorus sites, and 70 suspended sediment sites qualified for analysis. For the long-term trend interval, 44 total nitrogen sites, 18 total phosphorus sites, and 18 suspended sediment sites qualified for trend analysis.

The R package “EGRET” (ver. 3.0.7; Hirsch and De Cicco, 2015) contains the statistical model used for load and trend calculations called Weighted Regressions on Time, Discharge and Season (herein referred to as WRTDS; Hirsch and others, 2010). WRTDS is used to compute nitrogen, phosphorus, and suspended sediment daily concentration and load estimates. WRTDS was run using a dataset comprised of about 20 samples per year for each site and trend interval. After WRTDS was applied, a dynamic autocorrelation Kalman-filter (WRTDS-K) was used to adjust the estimated values based on the serial correlation of the residuals. WRTDS-K creates a daily time series of concentrations and loads that use the observed value when available and the serial-adjusted values otherwise. This approach provides more accurate concentration and load estimates, especially on sub-annual timescales (Zhang and Hirsch, 2019). These values were area-normalized to compute status estimates as yields.

The flow-normalized estimates from WRTDS were used to characterize trends over the two intervals. Estimates remove the influence of year-to-year variability of stationary streamflow on water quality and by doing so, provide a measure of temporal change in nitrogen, phosphorus, and suspended sediment loads independent of weather-driven streamflow variations. The flow-normalized values provide an indication of the effect of changing sources, delays associated with storage and transport of historical inputs, and any implemented management actions on stream water quality. Sites with flow-normalized loads that were lower at the end of the WY 1985–2020 or WY 2011–20 periods were classified as having “improved” conditions; whereas sites with flow-normalized loads higher at the end of each period were classified as having “degraded” conditions. WRTDS provides a likelihood estimate of a trend’s existence based on a probability score. A trend’s likelihood was categorized as “likely” (score was greater or equal to 0.67 and less than 0.9), “very likely” (greater or equal to 0.9 and less than 0.95), or “extremely likely” (greater or equal to 0.95). A site was classified as having no trend (meaning an improving or degrading trend is as likely to exist as it is not) if there is no discernable difference (likelihood estimate probability score between 0.33 and 0.67) between the flow-normalized loads in the start year and those in the end.

2.1.3. Results and Discussion

Status estimates for nitrogen yields ranged from 1.27 to 32.6 pounds per acre (lb/acre), total phosphorus yields ranged from 0.11 to 1.89 lb/acre, and suspended sediment yields ranged from 23.9 to 1,210 lb/acre (fig. 2; Mason and others, 2023). Higher total nitrogen yields were concentrated in the lower Susquehanna River Basin and Delmarva Peninsula, and lower yields were observed in the upper Susquehanna River and upper and lower Potomac River Basins and throughout the Rappahannock, York, and James River Basins. Similar high and low patterns for total phosphorus existed throughout the Susquehanna River Basin. Low yields were observed in the upper and lower Potomac River Basin, and in central Virginia saw a spread of high and low total phosphorus yields. Suspended sediment yields were highest towards the mouths of the Susquehanna and Potomac Rivers, and in central Virginia. The central Susquehanna River Basin, upper Potomac River Basin, and Eastern Shore recorded low yields relative to the rest of the Bay.

For short-term trends, 34 of 89 sites (38 percent) showed improvement for total nitrogen with yield reductions ranging from 0.06 to 2.26 lb/acre, and 37 sites (42 percent) had degrading trends with yield increases ranging from 0.02 to 2.62 lb/acre. Eighteen sites (20 percent) had no trend for total nitrogen. For total phosphorus, 31 of 70 sites (44 percent) had improving trends with yield reductions ranging from 0.007 to 0.31 lb/acre, and 16 (23 percent) had degrading trends with yield increases ranging from 0.009 to 0.82 lb/acre. Twenty-three sites (33 percent) had no trend for total phosphorus. For suspended sediment, 13 of 70 sites (18 percent) had improving trends and yield reductions ranging from 16.7 to 552 lb/acre; 32 sites (46 percent) had degrading trends and yield increases ranging from 9.11 to 2310 lb/acre. Twenty-five sites (36 percent) had no trend in suspended sediment (fig. 3).

Of the 44 sites analyzed, long-term trend results for total nitrogen show 24 sites improving, 15 sites degrading, and 5 sites with no discernable trend direction. The change in load ranged from −64.7 to 44.4 percent. Long-term trend results for the 18 sites analyzed for total phosphorus show 12 sites improving, 4 sites degrading, and 2 sites with no trend. The change in load ranged from −71.7 to 86.6 percent. Long-term trend results for the 18 sites analyzed for suspended sediment show 7 sites improving, 7 sites degrading, and 4 sites with no trend. The change in load ranged from −68.4 to 56.8 percent (fig. 3).

Additional information for each monitoring site is available through the USGS website “Water-Quality Loads and Trends at Nontidal Monitoring Stations in the Chesapeake Bay Watershed” (U.S. Geological Survey, 2016). This website provides State, Federal, local partners, and the public access to a wide range of data for nutrient and sediment conditions across the Chesapeake Bay watershed.

2.2.1. Data Compilation

A status analysis was conducted to describe the spatial variability of observed SC levels relative to background SC across the Chesapeake Bay watershed, and a trend analysis was conducted to quantify changes in SC over time. SC data for these analyses were obtained from a recently published SC inventory, which contains SC data compiled for the Chesapeake Bay watershed from 1901 to 2022 (Fanelli and others, 2023). These data originated from the WQP, which serves as a national repository for local, State, and Federal data collection efforts (U.S. Geological Survey and U.S. Environmental Protection Agency, 2021). The WQP also contains all available discrete data in NWIS.

Freshwater nontidal streams and rivers were the focus of these analyses, and sites with known tidal influence were removed from further analysis. This was done by evaluating and excluding sites that intersected areas near the Chesapeake Bay estuary that were designated by the National Oceanic and Atmospheric Administration as tidally influenced (Chesapeake Bay Program, 2004). Tidal influence was also manually verified by screening the remaining sites for high values of SC (1,000 µS/cm or greater at sites adjacent to the tidal zone), which may indicate tidal influence. The remaining sites were considered for subsequent status and trend analyses.

2.2.2. Analysis

All 278 sites identified in the tidal screening process were utilized for status and trend analysis. For the salinity analysis seasons were defined as December through February (winter), March through May (spring), June through August (summer), and September through November (fall). To accommodate these seasonal intervals and maximized the number of sites eligible for analysis across the watershed, a ‘year’ was defined as December 1 of the previous year to November 30 of the current year (for example the year 2008 was defined as December 1, 2007, through November 30, 2008). Hereafter within section 2.2, all years and year ranges refer to this defined time interval (for example, 2008–17 refers to December 1, 2007-November 30, 2017). At each site, SC status was computed as the median of the annual medians of 2015, 2016, and 2017. These 3-year median values were then compared to previously published background SC estimates that used SC observations from 2001–15 taken at reference sites and represent natural SC levels that would occur in the absence of anthropogenic inputs (Olson and Cormier, 2019). Background estimates were generated on a monthly timescale for years 2001–15 using a random forest model that predicted background SC using geological and climate variables across the United States for all National Hydrography Dataset (NHD) version 2.1 stream reaches (McKay and others, 2012). For the SC status analysis, the predicted monthly background SC estimates from Olson and Cormier (2019) were summarized as median annual background SC values for each of the years 2008–15, and a median of those values was computed to represent a long-term median annual background SC estimate. This long-term median annual background SC estimate was then compared to the observed 3-year median SC values.

For this comparison, each site was first associated to the National Hydrography Dataset version 2.1 flowlines to generate a list of stream segment unique identifiers (COMIDs) for each site. Next, the site COMID list was joined to the background SC dataset. Of the 278 sites, 269 (97 percent) had background SC data associated with their COMID; the remaining nine sites either did not have an associated COMID or there was no background SC reported in the 2019 study dataset for that COMID. Finally, the observed 3-year median was divided by the long-term median background SC to represent an “observed versus expected” metric. This metric characterizes the extent to which observed SC departed from background SC. SC departure classes were defined as “at or below the background SC estimate”, “slightly elevated” (1–2 times the background SC estimate), “moderately elevated” (2–3 times the background SC estimate), and “highly elevated” (greater than 3 times the background SC estimate).

Two trend analyses were conducted in R (ver. 4.2.3; R Core Team, 2023) using the compiled data: WRTDS (Hirsch and others, 2010) and Seasonal Mann-Kendall (SMK) trend tests (Helsel and others, 2020). WRTDS is a robust statistical model used to quantify water-quality trends, but it requires daily streamflow information, which limits its application to those sites with streamflow (discharge) data. SMK trend tests, however, only require water-quality information and can therefore be applied to a larger number of sites across the watershed. Results from both analyses can be found in Boyle and others (2025).

2.2.3. Results and Discussion

2.2.3.1. Status Results

The average long-term median background SC was 82 µS/cm for the 278-site network and ranged from 21 to 312 µS/cm (fig. 4). The lowest background SC estimates were in central Virginia and northern Pennsylvania. The highest background SC estimates were in carbonate regions in central Virginia and southeastern Pennsylvania. The average 3-year median for observed SC for all sites was 242 µS/cm across the 278-site network for the time period of December 1, 2014 to November 30, 2017 and ranged from 9 to 1,189 µS/cm across the watershed. The sites with the lowest 3-year median SC values were typically in central Virginia, in northwest Pennsylvania, and in the Delmarva Peninsula. The sites with the highest values were clustered in Maryland, especially throughout the Baltimore–Washington, D.C. metropolitan area, and in southeastern Pennsylvania. Departures from background SC follow a similar spatial pattern. Sites with values at or below background SC were mainly clustered in central Virginia (fig. 4). Only 43 sites (16 percent) had 3-year median SC values that were at or below background SC (table 3). Fifty-seven sites (21 percent) had slightly elevated values; these sites were clustered in central Virginia and western Pennsylvania. Sites that had moderately elevated or highly elevated values were throughout northern Virginia, Maryland, and Pennsylvania, and clustered in the Baltimore–Washington, D.C., metropolitan area, and in southeast Pennsylvania. These moderate and high departures from background SC were observed at 169 of the 269 sites with available background SC data (63 percent), indicating the effects of freshwater salinization are prevalent through much of the watershed. These results were consistent with a modeling study that estimated two-thirds of nontidal stream reaches were elevated above background SC in the Chesapeake Bay watershed (Fanelli and others, 2024).

Table 3.    

Counts and percentages of salinity sites with specific conductance (SC) status values for 2015–17 per SC departure class, and counts of salinity sites located where SC predicted background data are available.

[Data are from Boyle and others (2025). >, greater than; ≤, less than or equal to; NA, not applicable]

Table 3.    Counts and percentages of salinity sites with specific conductance (SC) status values for 2015–17 per SC departure class, and counts of salinity sites located where SC predicted background data are available.

SC departure class	Sites used in salinity status analysis
	Number	Percent
At or below background SC	43	16
>1 to ≤2 times background SC	57	21
>2 to ≤3 times background SC	73	27
>3 times background SC	96	36
No data available	9	NA
Total	278	100

2.2.3.2. WRTDS Trend Results

Ten-year SC trends from 2008 to 2017 were significant at most of the 35 sites in the WRTDS analysis; 21 sites (60 percent) had significantly increasing trends, and 6 sites (17 percent) had significantly decreasing trends (table 4). These results were consistent with other regional and national SC trend analyses that found general increases in SC (Kaushal and others, 2018; Bird and others, 2018; Baker and others, 2019). Trends were in the “highly likely” category for 14 of the 21 sites with significantly increasing trends, indicating higher confidence in those patterns. Sites with increasing trends were scattered throughout the Chesapeake Bay watershed (fig. 5). Sites with decreasing trends were largely found in the Potomac River Basin in Maryland and in central Virginia. In general, rates of change estimates by WRTDS were small; annual increases were below 5 µS/cm per year in all but one of the upward trending sites (fig. 6). The exception was a site in downtown Washington, D.C., at which the median annual SC was estimated to be increasing approximately 11 µS/cm per year during the trend period. Declines in SC at most of the six downward trending sites were also modest (less than −5 µS/cm per year). Annual seasonal change estimates indicated that greater changes in SC occurred in the winter and spring seasons at sites with increasing trends (fig. 7), which aligns with patterns observed in the Delaware River Basin (Rumsey and others, 2023). By contrast, the largest declines in annual seasonal change estimates occurred in the summer and fall seasons in sites with significant decreasing trends.

Table 4.    

Trend direction and significance categories for salinity sites across the Chesapeake Bay watershed analyzed using Weighted Regressions in Time, Discharge, and Season (WRTDS) and Seasonal Mann-Kendall (SMK) trend analyses, 2008–17.

[Data are from Boyle and others (2025). NA, not applicable]

Table 4.    Trend direction and significance categories for salinity sites across the Chesapeake Bay watershed analyzed using Weighted Regressions in Time, Discharge, and Season (WRTDS) and Seasonal Mann-Kendall (SMK) trend analyses, 2008–17.

Trend direction and significance	Sites analyzed with WRTDS		Sites analyzed with SMK
	Number	Percent	Number	Percent
Upward trend is highly likely	14	40	46	17
Upward trend is very likely	2	6	22	8
Upward trend is likely	5	14	24	9
No trend	8	23	152	55
Downward trend is likely	2	6	15	5
Downward trend is highly likely	4	11	13	5
Downward trend is very likely	0	NA	6	2
Total	35	100	278	100

2.2.3.3. Seasonal Mann-Kendall Trend Results

Significant trends were detected at slightly less than half of the SMK trend sites; increasing trends were detected at 92 sites (33 percent) and decreasing trends, 34 sites (12 percent; table 4). No trends were detected at the remaining 152 sites. Half of the detected increasing trends were in the “highly likely” category (46 sites), indicating high confidence in these trends. By contrast, only 38 percent of the detected decreasing trends were in the “highly likely” category (13 sites). Upward trending sites were scattered throughout the Chesapeake Bay watershed, following a similar pattern seen in the WRTDS results (fig. 5). There was also a cluster of upward trending sites with high trend certainty in the Baltimore–Washington, D.C., metropolitan area. Sites with decreasing trends, however, were only found in the southern portion of the watershed, typically south of the Pennsylvania border and especially in central Virginia (fig. 5).

Theil-Sen slope estimates provide annual estimates of rates of change in SC for those sites with significant trends. Most (61 sites or 66 percent) of the upward trending sites had small change estimates (increasing at a rate of less than 5 µS/cm per year), following the same pattern seen in the WRTDS results (fig. 6). However, 13 percent of the upward trending sites have change estimates greater than 10 µS/cm per year, suggesting rapidly worsening conditions at these sites, almost all of which were in the Baltimore–Washington, D.C., metropolitan area. Sites with moderate annual increases in SC (increases between 5 and 10 µS/cm per year) were also found scattered throughout Pennsylvania. By contrast, most (85 percent of the 34 sites) annual change estimates in sites with decreasing trends were small (less than −5 µS/cm per year), indicating modest declines over time.

The SMK and WRTDS trend results for the 35 sites that fit criteria for both analyses were compared. Results among the two methods generally agreed, and trend directions from the two methods were not found to be contradictory. More information on the comparison may be found in appendix 1.

In general, many of the increasing trends in SC co-occurred in areas with elevated SC (fig. 8). For example, 40 of the 93 SMK sites (43 percent) that have highly elevated 3-year median SC values also have a significantly increasing trend in SC (table 3). The WRTDS sites show a similar pattern: 8 of the 14 sites (57 percent) in the highly elevated SC departure class also have increasing SC trends. The highly elevated SC values at these sites result from 10 years of increasing SC. By contrast, almost all (91 percent) of the SMK sites that had 3-year median SC values at or below background SC also had either no change or decreasing SC trends, suggesting that low SC conditions have persisted at these sites. One possible reason for this is that these sites are in areas that are protected from land-use change or other anthropogenic disturbance. It is important to note that the distribution of trend sites among the departure classes is different between the two trend types; the SMK sites are more evenly distributed among the four departure classes, suggesting they cover a range of land-use settings. There are no WRTDS sites in the “at or below background” departure class, however, indicating little representation of undisturbed settings within those results.

The largest annual rates of increase and decrease were observed at sites with highly elevated SC (fig. 8), suggesting the conditions in these sites are highly dynamic. Rates of increase exceeded 20 µS/cm per year at a handful of sites, all of which had 3-year medians at or above 500 µS/cm. The same pattern is noted in sites with declining trends—the largest SC declines are also at sites with high 3-year median SC values. This could be due to point-source reduction efforts. Another trend study documented some decreasing SC trends in the region (Bowen and others, 2015), but the drivers of those patterns were not investigated. Additional studies could help fully understand the factors controlling these spatial and temporal patterns in SC.

2.3.1. Data Compilation

Continuous and discrete datasets can differ substantially in sample size, frequency, and spatial distribution, and these population dynamics are important considerations for determining status and trend approaches. Continuous stream temperature datasets with evenly spaced high-frequency time series (such as every 15 minutes) have become more common over the past decade. These data are often collected at fixed locations along a reach and are representative of diurnal and seasonal temperature fluctuations but are collected at a small number of sites that are often not spatially representative across the watershed. In comparison, discrete stream temperature datasets with unevenly spaced time series and infrequent measurements (6–8 per year) may be available over longer periods. However, discrete measurements may be less representative of site conditions, and diurnal and seasonal temperature fluctuations because measurements are not always collected at fixed locations or time intervals.

Data compilation and analyses were completed in R (ver. 4.2.1; R Core Team, 2022b). Stream temperature measurements were collated across the Chesapeake Bay watershed from NWIS, the WQP, and the USGS Aquarius Time-Series database (Aquarius) for consideration in status and trend analysis (fig. 9; Clune and others, 2023). Continuous unit values (high-frequency, evenly spaced time series) and daily aggregate data (minimum, maximum, and mean) of stream temperature were retrieved for USGS sites from NWIS using the R package “dataRetrieval” (ver. 2.7.11; Hirsch and De Cicco, 2015; U.S. Geological Survey, 2023). Depending on the site, daily stream temperature values can have long historical records. In contrast, many continuous unit values are only available after 2007. Stream temperature was measured with a sensor at fixed locations or depths. Sometimes, multiple overlapping time series are available if the sensor was relocated or replaced (Clune and others, 2023). These continuous unit values and aggregate daily data have been approved or considered provisional, revised, or estimated under USGS guidelines and procedures for measurement of stream temperature (Wagner and others, 2006). Any data point that was outside five standard deviations from the monthly mean was verified with the responsible USGS Water Science Center and represents infrequent or atypical stream temperature measurements that are deemed accurate.

Discrete values (infrequent, unevenly spaced time series) for stream temperature from multiple agencies stored in the WQP for streams within the Chesapeake Bay watershed were also obtained using the R package “dataRetrieval” (ver. 2.7.11; Hirsch and De Cicco, 2015). Additionally, discrete stream temperature values collected by the USGS while measuring streamflow were obtained from an internal pull of the Aquarius database. These stream temperature measurements were collected as an instrumentation check following techniques and standards for streamflow measurements (Turnipseed and Sauer, 2010). Measured accuracies are not available because these temperature measurements have not been approved under USGS guidelines and procedures for measurement of stream temperature (Wagner and others, 2006; U.S. Geological Survey, 2024). The measuring point of these discrete stream temperature measurements is not based on a fixed location or depth and may vary horizontally and vertically within the stream channel. Missing and unreasonable stream temperature discrete values outside the range of −5 and 40 °C were removed.

2.3.2. Analysis

Developing status and trends of stream temperature can be challenging, and there are advantages and limitations to using these data for analysis. Data bias and errors may be inherent, especially with multi-agency discrete data (Wilby and others, 2017). Temperatures can vary horizontally and vertically within the stream, and changing measurement points during collection could therefore influence results. Sites downstream of local dams, impervious areas, or water treatment outfalls can have altered thermal regimes and may present break points in the time series and introduce bias in regional models for stream temperature. Changes in the accuracy of instrumentation (mercury bulb thermometer to thermistor) during the period of record are often unknown. The collection times of discrete stream temperature data often exclude nights and weekends and may favor the morning or afternoon or may change over the period of record. Further, important metadata, such as sampling date and time, site description, and hydrologic condition, are often incomplete, inaccurate, or unavailable for discrete data. For example, 77 percent of Aquarius data were provided with quality control flags for stream temperature sample times that may be erroneous (outside working hours), not accurate (up to 4 hours off), or not available and only accurate to the day. Identifying such issues is not always feasible for multiple agency datasets like the WQP. Lastly, information on environmental conditions that can profoundly influence stream temperature, such as streamflow, is often missing from discrete stream temperature data.

Distinguishing whether a relation between two variables changes in a monotonic or nonmonotonic fashion can help determine the appropriate procedures and inferences for trend analysis (Wagner and others, 2017; Helsel and others, 2020). Monotonic trend approaches assume a gradual increase or decrease in stream temperature over time and have been used frequently for trend analysis of stream temperature (Ashizawa and Cole, 1994; Kaushal and others, 2010; Rice and Jastram, 2015). Trend analysis of annual or monthly mean stream temperature often employs parametric (linear regression) or non-parametric (Mann-Kendall with Theil-Sen slope) methods. Additional techniques that adjust for seasonality in monthly water-quality data (such as SMK) have been applied to stream temperature trend analysis.

These methods can be informative if the annual or monthly means are based on continuous data. Monthly and annual means of stream temperature computed from a small population of discrete data collected at varying points within the diurnal and seasonal cycle are seldom comparable to similar aggregate values based on routine, evenly spaced, continuous measurements over the same period. Nonmonotonic trend approaches, such as Seasonal-Trend decomposition using locally estimated scatterplot smoothing (LOESS) or WRTDS, have been developed to address the diurnal and seasonal patterns with environmental data that do not follow a unidirectional change over time (Cleveland and others, 1990; Hirsch and others, 2010). Although these methods have been widely used to detect trends in other water-quality data, they do not fully address the assumption of sample independence (serial correlation) inherent with continuous stream temperature time-series data, nor issues with the time of day or where in stream channel the measurement was made. Stream temperature trend approaches need to have the flexibility to capture nonlinear relationships and allow for inclusion of multiple predictors such as day, time of day, and streamflow while not assuming independence among measurements. Given the characteristics of the data and this consideration, different approaches were used to estimate the status and trends for the continuous and discrete temperature datasets.

2.3.2.1. Continuous Data

Generalized additive models (GAMs) were the chosen trend analysis method for this study because they are ideal for capturing a smooth, nonmonotonic trend estimation and nonlinear effects of covariates for datasets that have a limited number of sites with extended continuous time series (Murphy and others, 2019; Yang and Moyer, 2020). Continuous stream temperature records from NWIS datasets (daily mean values) published by Clune and others (2023) were used for the GAMs (fig. 9). The sites selected for this study had stream temperature collection data from at least one of four contemporary periods (WY 2013–20, WY 2013–21, WY 2013–22, and WY 2014–22) during which at least 300 daily values per year were available. Though differing trend start dates prevent direct trend comparison (Helsel and others, 2020), the use of these multiple contemporary periods maximizes the site count, spatial coverage, period of record and sample size across the Chesapeake Bay watershed. Using longer trend intervals (those greater than 10 years) or excluding contemporary periods that may be missing only 1 year at the beginning or end of the record would significantly reduce the spatial extent of sites across the watershed. Datasets with less than 8 years of data are considered too short for separating trends from seasonal fluctuations in water quality (Murphy and others, 2019). Thirty-one sites met these criteria with an average of 346 daily value measurements per site (fig. 10). Drainage areas ranged from 1.01 to 11,220 square miles and had a median size of 49.1 square miles. No duplicate or erroneous values were apparent. Continuous data are considered to be representative of mean fluctuations in stream temperature over a 24-hour period.

Continuous stream temperature trends were fit for each selected site using a GAM that uses a penalized regression-spline smoothing approach (Yang and Moyer, 2020). The response of continuous stream temperature (daily values) by site was expressed as a smooth function of covariates in the following form:

$g\left(T_{continuous}\right)=B_o+s\left(\log (Q\right))+s\left(year\right)+s\left(doy\right)+t_i\left(\log (Q\right),doy)$

(1)

where

g(T_continuous)

is stream temperature (daily values),

B_o

is the modeled intercept,

s(log(Q))

is the smooth function by log-transformed streamflow,

s(year)

is the smooth function by time (decimal year),

s(doy)

is the smooth function by season (day of year), and

t_i(log(Q), doy)

is the tensor product interaction of covariates of log-transformed streamflow and day of year (doy).

The R package “mgcv” was utilized to fit the GAMs for each site (ver. 1.8-42; Wood, 2017). Based on the effective degrees of freedom for each response variable, a k-value (number of knots) of 10 was used for the log-transformed streamflow and time smooth functions with cubic spline, and a k of 20 was used for the seasonal smooth function (Yang and Moyer, 2020). The direction and magnitude of the trend was calculated as the difference between stream temperature from the start and end of the trend interval. Trend results at the 95-percent confidence level were computed for each site and likelihood values were calculated to provide a nuanced inference and gradual estimation of probabilities (p-values; McBride, 2019). The upward or downward trend likelihood categories included “unlikely” (less than 33 percent), “about as likely as not” (33–67 percent), “likely” (67–90 percent), “very likely” (90–95 percent), and “extremely likely” (95–100 percent; McBride, 2019). The performance of each GAM by site was based on adjusted coefficients of determination ($R_{adj}^2$). Significance of covariate effects on water temperature were based on the p-value provided by the Wald statistics in the “mgcv” package (ver. 1.8-42; Wood, 2017). Additionally, the effective degrees of freedom (EDF) provided by the GAMs was used to quantify the degree of nonlinearity of the smooth functions.

In addition to computing trends in continuous stream temperature, these sites and daily values were used to provide the current status of stream temperature across the Chesapeake Bay watershed. The status of stream temperature was computed for each site as the departure of 2022 temperatures, in degrees Celsius above or below the mean annual temperature for the reference period (trend interval). The sparse distribution of continuous sites selected was too heterogenous in topography, land use, and ecoregion to develop an effective, regional status based on a specific species benchmark (such as Salvelinus fontinalis [Mitchill, 1814; brook trout]), disparate reference conditions, or contrasting statewide water-quality standards and criteria (Hawkins and others, 2010). Results from the status and trend analysis of continuous stream temperature data can be found in Boyle and others (2025).

2.3.2.2. Discrete Data

Datasets that have an abundant number of sites with infrequent and limited collection of discrete data at each site were evaluated with linear mixed models (LMMs) to produce a single, watershed-wide trend that accounts for seasonal and diurnal effects (Bolker and others, 2009). Discrete stream temperature records from the WQP and Aquarius datasets published by Clune and others (2023) were used in LMMs (fig. 9). Stream water temperature sites with discrete data were selected to maximize the spatial coverage across the Chesapeake Bay watershed and to provide the longest possible record and comparability in sample size for the trend period. The selected sites had data collected from calendar years 1985 to 2022 and had no more than 5 years of missing data. Within the Aquarius database, data before 1985 within the Aquarius database lacks important metadata (such as the time of day a sample was collected) needed for analysis. Additionally, data during the 1985–2022 trend interval provides good spatial representation and consistent sample size primarily because of the inception of the Chesapeake Bay Water Quality Monitoring program in 1985. The selected discrete stream temperature records included 116,535 measurements from 369 sites (fig. 11). On average, the number of stream temperature observations ranged from approximately 3 to 11 per year. Manual geospatial inspection of site locations revealed duplicate values for 49 sites, which were excluded from subsequent analysis. An additional 608 measurements were excluded because reported stream temperatures were below 0 °C. The final dataset for analysis of discrete data included 94,148 measurements across 320 sites with an average of 294 measurements per site.

Discrete stream temperature data ranged from 0.1 to 39.0 °C, and more than 90 percent of the samples were between 1 and 25 °C (fig. 11A). The dataset equally encompassed all seasons, as seen through the even distribution of samples across all days of the calendar year (day of year; fig. 11B) almost exclusively represented daytime hours of collection between 9 a.m. and 5 p.m. (eastern time; fig. 11C) and included samples each year from 1985 to 2022 (fig. 11D). The dataset contained an average of 2,478 measurements per year. The expected pattern of variation across day of year in discrete stream temperature data was apparent (fig. 12); the highest values were observed during summer months and lowest values, during winter months. This nonlinear association between temperature and day of year justified the incorporation of a quadratic term for day of year in the LMMs.

Using the screened discrete data, LMMs with fixed and random effects were used to characterize a watershed-wide trend in stream water temperature. Fixed-effect covariates included day of year, day-of-year squared to account for nonlinear seasonal effects, hour of day to account for diurnal variation, and year to assess temporal trends. Covariates were scaled to a mean of 0 and standard deviation of 1 (z-score transformation) to facilitate comparison of covariate effects. The model, which includes a random effect (intercept) by site to account for unmeasured effects of elevation, land use, stream volume, and other local effects on baseline water temperatures, was in the following form:

$f\left(T_{discrete}\right)=z\left(doy\right)+z({\text{doy}}^2)+z\left(year\right)+z\left(hour\right)+RE\left(site\right)$

(2)

where

f(T_discrete)

is discrete stream temperature,

z(doy)

is the season variation by day of year (z-score transformation),

z(doy²)

is the season variation by day of year (z-score transformation) squared,

z(year)

is the annual variation by decimal year (z-score transformation),

z(hour)

is the diurnal variation by hour of day (z-score transformation), and

RE(site)

is the location specific random effect (intercept) by site.

Functions in the R package “lme4” (ver. 1.1-32; Bates and others, 2015) were used to fit LMMs, and functions in R package “MuMIn” (ver. 1.47.5; Bartoń, 2010) were used to evaluate all pairwise additive combinations of candidate covariates based on Akaike information criterion (AIC). Interpretation of the goodness of fit of the best model was based on R² values for the conditional (fixed effects only) and marginal (fixed and random effects) coefficients of determination following Nakagawa and Schielzeth (2013) with functions in R package “performance” (ver. 0.10.3; Lüdecke and others, 2021).

Interpretation of the significance of covariate effects on stream temperature, including the time trend represented by z(year) in the model (eq. 2), was based on the departure of 95-percent confidence intervals from zero. Functions in the R package “glmmTMB” (ver. 1.1.7; Brooks and others, 2017) were used to estimate covariate confidence intervals, and functions in R package “sjPlot” (ver. 2.8.15; Lüdecke, 2013) were used to visualize the variation in random effects (random intercepts) among sites. Results from the analysis of discrete stream temperature data can be found in Boyle and others (2025).

2.3.3. Results and Discussion

In this section, the status of stream temperature is presented as the departure of 2022 water temperature from the mean of annual temperatures for the trend interval (WY 2013–20, WY 2013–21, WY 2013–22, or WY 2014–22, depending on site). Next, the site-specific trend results derived from the continuous stream temperature data and the generalized additive models are summarized. Finally, the watershed-wide temperature trend derived from discrete stream temperature data and the LMM is presented.

2.3.3.1. Status

WY 2022 was the second warmest year for stream temperature on average for sites with a 9- or 10-year reference period of record that starts in WY 2013 or WY 2014 and ends in WY 2022 (24 of 31 sites; table 5). The average mean annual temperature for the select sites in WY 2022 ranged from 9.66 to 16.14 °C. The departures were positive at all sites and ranged from 0.01 to 0.85 °C above the average for the period. These departures from the average for the period for each site differed spatially and the seasonal fluctuation in stream temperature represents a substantial fraction of the annual variability (table 5, fig. 13).

2.3.3.2. Continuous Data Trends

GAM methodologies that use the covariates of streamflow, season, year, and interaction term (streamflow and season) were selected to develop trends across sites with daily mean values derived from continuous stream temperature in the Chesapeake Bay watershed. Overall, the GAMs performed well across all 31 continuous data sites. The $R_{adj}^2$ ranged from 0.86 to 0.96. Nearly all smooth functions for covariates except two were significant in the fitted GAMs. Additionally, the effective degrees of freedom were consistent, providing a balance between model complexity and fit (table 5).

Table 5.    

Generalized additive model (GAM) comparison of 31 sites with continuous (daily) stream temperature in the Chesapeake Bay watershed.

[Data are from Boyle and others (2025). Values are not significant unless otherwise marked by asterisks: *, value has a significance level of 0.01; and **, value has a significance level of 0.05. USGS, U.S. Geological Survey; $R_{adj}^2$, adjusted coefficient of determination; Geo mean, geometric mean; s(log(Q)), streamflow (log); s(doy), season (day of year); s(year), time period (decimal year); t_i(log(Q),doy), interaction term for streamflow and season; °C, degrees Celsius; WY, water year; NA, not applicable]

Table 5.    Generalized additive model (GAM) comparison of 31 sites with continuous (daily) stream temperature in the Chesapeake Bay watershed.

USGS station number	Model		Effective degrees of freedom				Trend						Status
	$R_{adj}^2$	Geo mean1	s(log(Q))	s(doy)	s(year)	ti(log(Q),doy)		Trend interval2	Trend direction and likelihood3	Likelihood value3	Difference of annual mean from the start to the end of the trend interval, in °C	Percent change in temperature per year	Annual mean of trend interval, in °C	Departure of 2022 annual mean, in °C above or below the average annual mean of 2013–22 or 2014–22 trend interval
01540500	0.96	13.14	7.98**	14.09**	6.83**	7.93**		WY 2013–22	Upward trend is likely	0.79	0.34	0.26	12.67	0.85
01542500	0.95	12.38	3.59**	11.11**	2.66**	2.72**		WY 2013–22	Downward trend is extremely likely	0.975	−0.89	−0.72	12.30	0.41
01544500	0.94	10.27	5.32*	11.36**	7.00**	9.27**		WY 2013–22	Upward trend is about as likely as not	0.525	0.02	0.02	10.22	0.81
01545600	0.94	9.49	2.82**	11.33**	7.05**	8.74**		WY 2013–22	Upward trend is extremely likely	0.975	0.59	0.62	9.45	0.62
01549700	0.95	12.01	1.00**	11.10**	6.30**	8.69**		WY 2013–22	Upward trend is about as likely as not	0.66	0.18	0.15	11.85	0.70
01573695	0.93	12.39	1.00**	10.98**	1.00*	8.08**		WY 2013–22	Upward trend is extremely likely	0.975	0.7	0.56	12.36	0.31
01573710	0.93	12.72	4.68**	11.05**	1.00*	6.53**		WY 2013–22	Upward trend is extremely likely	0.975	0.75	0.59	12.68	0.30
015765195	0.88	12.12	4.16**	10.20**	1.00	7.50**		WY 2013–22	Upward trend is likely	0.785	0.26	0.22	12.25	0.18
01581920	0.86	9.68	7.26**	13.42**	1.00**	11.51**		WY 2013–22	Upward trend is extremely likely	0.975	0.84	0.87	9.66	0.13
01595800	0.94	10.19	6.23**	13.89**	3.99**	11.50**		WY 2013–22	Upward trend is likely	0.72	0.38	0.37	10.36	0.01
01596500	0.92	9.72	1.00**	10.80**	2.23**	4.61**		WY 2013–22	Upward trend is extremely likely	0.975	1.09	1.12	9.99	0.79
01597500	0.94	9.71	7.80**	14.15**	1.89	11.44**		WY 2013–22	Upward trend is very likely	0.95	0.94	0.97	9.84	0.23
01645704	0.92	13.69	4.42**	11.22**	5.92**	7.36**		WY 2013–22	Upward trend is likely	0.79	0.31	0.23	13.67	0.35
01645762	0.92	13.3	2.64**	11.14**	6.32**	6.95**		WY 2013–22	Upward trend is likely	0.725	0.22	0.17	13.25	0.36
01646305	0.91	13.6	2.74**	10.76**	6.13**	7.24**		WY 2013–22	Upward trend is likely	0.835	0.39	0.29	13.56	0.42
01648010	0.93	13.89	1.00**	10.20**	1.58**	7.21**		WY 2013–22	Upward trend is extremely likely	0.975	1.19	0.86	14.18	0.34
01649190	0.92	13.21	3.51**	11.33**	5.97**	5.40**		WY 2013–22	Upward trend is likely	0.675	0.19	0.14	13.17	0.43
01649500	0.93	14.89	4.02**	10.89**	5.68**	7.06**		WY 2013–22	Upward trend is about as likely as not	0.59	0.1	0.07	14.85	0.15
01650800	0.91	13.94	3.13**	10.50**	1.00**	6.60**		WY 2013–22	Upward trend is extremely likely	0.975	0.85	0.61	13.96	0.36
01656903	0.91	14.82	4.23**	10.98**	4.96**	8.20**		WY 2013–22	Upward trend is extremely likely	0.975	1.56	1.17	14.88	0.41
01581752	0.91	12.63	4.02**	10.73**	3.12**	6.91**		WY 2014–22	Upward trend is extremely likely	0.975	1.87	1.65	12.84	0.05
01651770	0.95	15.72	5.80**	11.59**	1.21*	8.83**		WY 2014–22	Upward trend is extremely likely	0.975	0.87	0.61	15.79	0.35
01651800	0.92	15.04	3.21**	10.64**	1.71*	8.30**		WY 2014–22	Upward trend is extremely likely	0.975	0.89	0.59	15.19	0.08
01654500	0.92	13.78	1.49**	11.43**	2.66**	7.53**		WY 2014–22	Upward trend is extremely likely	0.975	1.23	0.99	14.13	0.44
02011500	0.94	13.18	1.00	11.24**	4.72**	6.87**		WY 2013–21	Downward trend is likely	0.75	−0.43	−0.33	NA	NA
02042500	0.91	16.67	1.00**	10.79**	2.24**	7.56**		WY 2013–21	Upward trend is very likely	0.925	0.85	0.51	NA	NA
02011400	0.93	13.3	4.33**	10.98**	5.11**	4.75**		WY 2013–21	Upward trend is extremely likely	0.975	1.05	0.88	NA	NA
01608500	0.94	14.67	4.11**	10.71**	5.36**	8.10**		WY 2013–20	Upward trend is about as likely as not	0.65	0.18	0.15	NA	NA
01610400	0.92	11.41	1.31**	10.56**	5.86**	2.46**		WY 2013–20	Upward trend is about as likely as not	0.56	0.08	0.08	NA	NA
01611500	0.95	14.47	3.67**	11.24**	5.47**	7.40**		WY 2013–20	Upward trend is likely	0.78	0.42	0.36	NA	NA
01646000	0.92	13.4	2.81**	9.71**	4.31**	7.30**		WY 2013–20	Upward trend is likely	0.8	0.53	0.49	NA	NA

¹

The central tendency of change.

²

The period of record used for trend analysis.

³

The probability description and value that indicates increasing or decreasing trends.

The semidecadal (8–10 years) trend estimates in stream temperature were calculated for four contemporary periods (WY 2013–20, WY 2013–21, WY 2013–22, and WY 2014–22; table 5). Increasing trends in stream temperature were “likely” to “extremely likely” for 79 percent of the sites across the Chesapeake Bay watershed, and only two sites indicated downward trends (fig. 14). For sites with increasing trends, total warming from the beginning to the end of the trend interval ranged from 0.19 to 1.09 °C (table 5).

GAMs provide the ability to further examine temporal change and the effect of each covariate on the trend. For example, there is an extremely likely and significant increasing trend of 1.56 °C at Flatlick Branch above Frog Branch at Chantilly, Virginia (USGS 01656903) from 2012 to 2022 (fig. 15A). The most significant rate of change occurred at this site from 2014 to 2017 (fig. 15B). Besides the temperature variation during the year (approximately 20 °C), rising stream temperature at this site is only slightly affected by streamflow and more influenced by the time (year) smooth function (figs. 15C–E).

These trend results for continuous stream temperature (table 5; fig. 14) demonstrate that the use of the novel GAMs approach by Yang and Moyer (2020) is suitable for wider use in the Chesapeake Bay watershed. Contemporary studies that also use continuous daily mean data have shown similar increasing trends in stream temperature for rivers nationally (Tassone and others, 2022; Zhi and others, 2023). Continued warming of streams over extended periods may surpass the thermal limits of aquatic species, especially those in cold water streams, and is a concern for native brook trout and other fisheries. Additionally, there is growing concern the implementation of best management practices for specific water-quality goals (for example, nutrient and sediment reduction) may be having adverse impacts on stream temperature and be acting as “heaters” rather than “coolers” to streams (Batiuk and others, 2023). Future work could utilize continuous daily mean, minimum and maximum values from the USGS and other agencies to better evaluate the full range of thermal regime (magnitude, variability, frequency, duration, and timing) that impact warm and cold water biota (Arismendi and others, 2013). Ultimately, the use of covariate effects in a GAM can be coupled with other analyses of climate, land use, and additional factors to better detect and explain the spatial and temporal drivers of change in stream temperature throughout Chesapeake Bay watershed.

2.3.3.3. Discrete Data Trends

LMMs using discrete data were utilized to produce a watershed-wide trend that accounts for seasonal and diurnal effects. Comparison of alternative combinations of covariates revealed one LMM for stream temperature that strongly outperformed all others (table 6). This model included an effect of hour of day, day of year, day-of-year squared, year, and random intercepts by site. The importance of the year effect was evident because the best-performing model (M1) included all covariates, whereas the next-best model (M2) lacked a year effect and had a much higher AIC value than M1, indicating M1 had a better fit. M1 accounted for 73.0 percent of the observed variation in stream temperature by fixed effects alone (marginal R²) and 77.2 percent of the measurement variation with the inclusion of fixed and random effects (conditional R²).

Table 6.    

The four top-performing linear mixed models of stream temperature in the Chesapeake Bay watershed based on all pairwise additive combination of covariates.

[Data are from Boyle and others (2025). All models included random intercepts for each site and a global y-intercept of 13.50. Coefficients are given for model covariates with degrees of freedom (df) and Akaike information criterion (AIC) values for candidate models. ΔAIC; difference in AIC value between each model and the top model (M1), NA; not applicable]

Table 6.    The four top-performing linear mixed models of stream temperature in the Chesapeake Bay watershed based on all pairwise additive combination of covariates.

Model rank	Covariate				df	AIC	ΔAIC
	Hour	Day of year	Day-of-year squared	Year
M1	0.45	26.14	−24.79	0.14	7	511,581	0.00
M2	0.45	26.14	−24.79	NA	6	511,705	123.24
M3	NA	25.97	−24.63	0.12	6	512,819	1,237.38
M4	NA	25.97	−24.63	NA	5	512,917	1,335.49

Inspection of confidence intervals for covariates in M1 revealed a significant warming trend in stream temperature from 1985 to 2022 (table 7). Specifically, the 95-percent confidence interval for the model coefficient on the year effect (0.14) ranged from 0.11 to 0.16; because these values do not encompass zero, they indicate a significant warming trend over the years and sites included in the discrete stream temperature dataset. The magnitude of seasonal effects exceeded the magnitude of hour effects or year effects in M1 (greater absolute value of scaled coefficients; tables 6, 7). Estimated increases in mean annual stream temperatures across all sites were within 1 °C (fig. 16).

Table 7.    

The 95-percent confidence intervals defining the effect direction and magnitude for covariates in M1, the top-ranked model for stream temperature in the Chesapeake Bay watershed.

[Data are from Boyle and others (2025), Model rankings shown in table 6. Covariates were scaled to facilitate comparison]

Table 7.    The 95-percent confidence intervals defining the effect direction and magnitude for covariates in M1, the top-ranked model for stream temperature in the Chesapeake Bay watershed.

Model parameter	2.5-percent level	97.5-percent level
Global y-intercept	13.32	13.67
Day of year	26.05	26.24
Day-of-year squared	−24.89	−24.70
Hour	0.43	0.48
Year	0.11	0.16

Prior research that also used discrete data reported increasing stream temperatures in the headwaters of the Chesapeake Bay watershed (Kaushal and others, 2010; Rice and Jastram, 2015). Results from this study support previous findings that stream temperatures in the Chesapeake Bay watershed are increasing and demonstrate the utility of discrete data for temporal trend analysis of water temperature. Additional research could help understand spatial variation in warming trends within the study area, including analysis of riparian cover, land-use change, and effects of geology on groundwater-stream water interactions that may regulate stream sensitivity to atmospheric change (Snyder and others, 2015; Hitt and others, 2023; Kessler and others, 2023).

The ecological implications of rising stream temperatures are a concern to resource managers throughout the Chesapeake Bay watershed. WY 2022 was the second warmest year for stream temperature on average for continuous sites with a 9- or 10-year period of record (WY 2013–22 or WY 2014–22). Increasing trends (0.19–1.09 °C) in stream temperature were “likely” to “extremely likely” for 79 percent of the continuous sites across the Chesapeake Bay watershed, and only two sites indicated downward trends. Additionally, estimated increases in mean annual water temperatures across all discrete sites were within 1 °C. Identifying the magnitude and geographic variability of warming stream temperatures across the watershed is the first step in identifying problematic areas and supporting healthy aquatic ecosystems into the future. Recent advances in the public availability of discrete monitoring data and increased use of continuous monitoring data offer opportunities to better assess thermal regimes in freshwater systems. These results will provide a foundation for more routine and systemic tracking of status and trends of stream temperature for use by environmental managers in the Chesapeake Bay watershed.

2.4.1. Data Compilation

In June 2019, the USGS contacted Federal and State agencies across six States (Delaware, Pennsylvania, Maryland, New York, Virginia, and West Virginia) and the District of Columbia (Washington, D.C.) to request available data on toxic contaminants within the Chesapeake Bay watershed. In addition, NWIS and the WQP were queried to retrieve relevant toxic contaminant data. The purpose was to create an inventory of available data across multiple jurisdictions within the Chesapeake Bay watershed. Data compiled for the inventory were limited to only sites where Hg, PCBs, or organochlorine pesticides were collected with sufficient metadata such as sample media, analytical method, sample date, sampling coordinates, and frequency of collection. These toxic contaminants were selected because of their toxicity, prevalence, and TMDL enforcement within the watershed. The related analytical results (concentrations) were not retrieved or used in this analysis because the focus was on data availability.

Data were formatted into two consolidated datasets in a USGS data release (Banks and others, 2022). The “detailed” dataset contains a record of each sample by unique site number, constituent (PCB, pesticide, and Hg), and sampling date. If available, each record includes the site name, data source (provider of the data), location (latitude and longitude), parameter code, analytical method description, collecting agency, and media type. It does not include the analytical results. The other dataset summarizes basic information associated with each unique site: the sampling start and end date, an indication of whether the location was sampled multiple times, the site’s location, the data source, the collecting agency, and the number of sampling records for PCBs, pesticides, and Hg (Banks and others, 2022). Additional information on the agencies and consolidation methods applied to create the “detailed” dataset are provided in the metadata for the data release. When available, total PCB concentrations were also compiled from NWIS, the WQP, and the data provided by Federal and State agencies across Chesapeake Bay watershed. These data were consolidated in a separate table to provide a tool to access comparable PCB data within the Chesapeake Bay watershed. A detailed discussion on compatibility of PCB data and why only total PCBs were compiled can be found in appendix 2 of this report.

2.4.2. Analysis

Spatial distribution of sites, number of records, sample frequency, and the period of record for each site provided the basis for assessing the toxic contaminant inventory (Banks and others, 2022). Analytical method compatibility for PCBs between sampling periods at a given location was also evaluated.

Most sites had few samples collected or samples collected at irregular intervals insufficient for trend analysis. An evaluation was performed to locate sites with repeat sampling events where trend analysis could occur if sampling continued. To identify sites that may qualify for potential future trend analysis, a decision tree was developed based on the available data in the contaminant inventory. The evaluation was limited to PCBs because of data availability and to avoid duplicating work completed by Willacker and others (2020) for Hg in the Chesapeake Bay. Pesticides were not included in this evaluation because of the non-uniformity in pesticides analyzed across the watershed and a priority on pollutants with a TMDL. Future assessments may be possible to assess status and trends for pesticides by evaluating the number of toxic pollutants in a given watershed over time.

The following criteria were used to identify sites where PCB trend analysis could occur in the future if data collection efforts continue:

2.4.3. Results and Discussion

The inventory created for this effort contains information about data availability for pesticides, PCBs, and Hg, and about concentration data for PCBs, from 1938 through 2019 for various sites across the Chesapeake Bay watershed. A comparison of the number of records compiled by media type (table 8) found that more records for pesticides and PCBs are available in NWIS than in State data sources. Some States, such as Pennsylvania, utilized the USGS for monitoring studies, which resulted in most records being in NWIS instead of separate State databases. The number of records available associated with a single sampling point varies by reporting agency and data source. For example, each individual pesticide and media in NWIS has a separate parameter code generating a unique record. A State agency record may only contain total PCBs for a given sample generating one record.

Different distributions of data were observed when comparing PCB, Hg, and pesticides by sample media and data source. Most PCB records were available from States (12,546 records), of which 71 percent were biological samples, 21 percent were sediment samples, and 7 percent were water samples. This is consistent with EPA guidance for the development of PCB TMDLs that utilize biological samples (fish tissue) to assess the potential impacts to human health, and derive empirical bioaccumulation factors from water, sediment, and fish concentrations measured within the impacted watershed (U.S. Environmental Protection Agency, 2011). State records for Hg also contained more biological data (79 percent) than any of the other media. Records of Hg water samples were only reported within NWIS (U.S. Geological Survey, 2023). Pesticide data were primarily from NWIS and consisted mostly of water samples.

Sampling by State and local jurisdictions for toxic contaminants is primarily focused on assessing the risk to human health and therefore targets the media responsible for the exposure pathway for the specific contaminant. These exposure pathways may not fully encompass the potential ecological impacts for these contaminants or bioaccumulation models, making it difficult to assess broader impacts (Fuchsman and others, 2006; Rattner and McGowan, 2007; Gobas and Arnot, 2010). For example, EPA guidance for fish sample preparation differs by fish species to reflect how that species is typically consumed by people (U.S. Environmental Protection Agency, 2000).

A total of 89 sites met criteria for potential future trend analysis for PCBs (table 9; fig. 17). The most common media at these sites was biological, but some water samples were taken in the Potomac River and James River Basins. Future trend analysis could consist of evaluating PCB concentrations in fish tissue collected within the major drainage basins identified in table 9. Where applicable, additional factors such as land use could be incorporated to assess significant differences in PCB recovery of the watersheds. Evaluating statistical significance will be difficult because of the low number of records available at each site. Aggregating samples by watershed could increase the sample size but could increase variability in levels. PCB contamination and recovery is very site specific even within a small watershed. For example, a 5-acre lake on the St. James River was treated with activated carbon to reduce PCB bioaccumulation in fish (Patmont and others, 2020). This remedial action was successful; however, the reduction in fish tissue concentration was only observed in resident fish within the lake. Similarly, elevated levels of PCB contamination can be isolated to specific tributaries within a watershed. Beaverdam Creek (Washington, D.C.) within the Anacostia River watershed is one of the five primary tributaries and responsible for most of PCB inputs to the main stem of the river (Wilson, 2020; Lombard and others, 2023).

2.5.1. Data Compilation

All data compilation and analyses were completed in R (ver. 4.4.0; R Core Team 2024). Streamflow daily value data were acquired from NWIS (U.S. Geological Survey, 2023) using the R package “dataRetrieval” (ver. 2.7.16; Hirsch and De Cicco, 2015) for 587 streamgages in the Chesapeake Bay watershed that have records from WY 1986 through the close of WY 2022. Within this interval, site selection criteria required 37 complete water years of streamflow data to ensure accurate estimation of streamflow trends, resulting in a final selection of 211 qualifying sites.¹

2.5.2. Analysis

To analyze streamflow status and trends, six metrics of streamflow extrema describing high- and low-flow magnitude, frequency, and duration were calculated for each year in the 37-year interval (table 10). These metrics provide a representation of distinct flow regime components and have been linked to stream biological condition (Carlisle and others, 2017; Eng and others, 2017, 201952). Though all sites had 37 complete years of data, some sites had fewer than 37 years of low- or high-flow duration values because some years did not have any daily values exceeding the low- or high-flow percentile threshold. To calculate a streamflow status for each metric, the percentage difference between a given metric’s value in WY 2022 and the same metric’s WY 1986–2022 average value was used to calculate streamflow status for each metric. Trends in streamflow metrics were estimated using GAMs for each site with metric values available for at least 30 years in the 37-year interval, because the nonlinear qualities of GAM responses support characterization of changes in extrema over time. GAMs were fit using the R package “mgcv” (ver. 1.9-1; Wood, 2017) with a smoothed year covariate and compared to a null model with no year covariate to determine trend significance. A trend was considered strong and significant if the year covariate p-value was less than 0.05 and the AIC value was more than 7 points smaller than found in the null model. Conversely, failure to meet these criteria were interpreted as no trend. Results from the status and trend analysis of streamflow data can be found in Boyle and others (2025).

2.5.3. Results and Discussion

Low- and high-flow duration status values could not be calculated for 46 sites; 43 sites did not have any low-flow events (flows below the 10th percentile of all daily value site flows in WY 1986–2022), and 3 sites did not have any high-flow events (flows above the 90th percentile of all daily value site flows in WY 1986–2022) in WY 2022. For most sites, high-flow frequency was above average and high-flow magnitude and duration were below average, indicating that though high flows were more common in WY 2022, those high flows were shorter in duration and generally smaller in magnitude than the 37-year average (fig. 18). Most sites had above average low-flow magnitude and below average low-flow duration and frequency, indicating that low flows in WY 2022 were not as low as the 37-year average and occurred less often and over fewer total days than the 37-year average of low flows (fig. 18). During WY 2022, above average high-flow frequencies and below average high-flow magnitudes and durations were most prevalent in the lower portion of the watershed. Most above average low-flow frequencies were associated with sites in the lower portion of the watershed, as were most sites with above average low-flow magnitudes. Most sites where below average low-flow durations were observed, however, were in the upper portion of the watershed (fig. 18).

For most sites across the watershed, values of high- and low-flow metrics were extremely variable and had no strong directional trends over time. Low-flow metrics had a greater number of significant trends overall than high-flow metrics. Trends in low-flow metrics revealed increases in low-flow magnitude, decreases in low-flow frequency, and few trends in low-flow duration. Of the few sites with significant trends in high-flow metrics, most had increasing trends in high-flow magnitude and duration (table 11; fig. 19). Sites with significant trends in streamflow metrics were sparse and scattered across the watershed. However, sites with increasing trends in low-flow magnitude and decreasing trends in low-flow frequency were clearly clustered in the Baltimore-Washington, D.C., metropolitan area, especially in tributaries of the Potomac River (fig. 19), indicating that low-flow events in this area have become less frequent and less extreme in magnitude over time.

Natural fluctuations in flow are important for maintaining healthy stream biological assemblages (Poff and others, 1997). Climate change and anthropogenic activities such as urban development or water withdrawals can drive unnatural variations in streamflow that harm aquatic life (Döll and Zhang, 2010; Brown and others, 2013). Extreme low flows can cause declines in sensitive macroinvertebrate taxa and taxa richness and decreases in fish growth rates, survival, and abundance—both directly through reductions in available swift water habitats and, more often, indirectly through increases in water temperature and SC during low-flow conditions (Dewson and others, 2007; Miller and others, 2007; Walters, 2016). Extreme high flows also influence stream biological assemblages by physically washing organisms downstream, scouring away algal food sources, destroying habitat, and introducing pulses of sediment and nutrients, all of which lead to decreasing macroinvertebrate densities, richness, and abundance and altered fish population densities and structure (Scrimgeour and Winterbourn, 1989; Carline and McCullough, 2003; Chattopadhyay and others, 2021). At most sites across the watershed, values of high- and low-flow metrics were extremely variable with no strong directional trends over time. The lack of general watershed-wide patterns in streamflow metrics (assuming these are not the result of Type II errors [when analysis fails to detect an existing trend]) suggests that low or high-flow metrics at the few sites with significant change may be driven by site-specific factors rather than watershed-wide factors such as climate. Future investigation of site-specific drivers, incorporation of climatic variables into the modeling framework, and evaluations of additional streamflow metrics will aid in a comprehensive understanding and maintenance of freshwater flows in the Chesapeake Bay watershed.

2.6.1. Data Compilation

Two distinct but complementary analyses using different data sources were undertaken to evaluate the status and trends in hydromorphology across the Chesapeake Bay watershed: (1) status and trends in multi-jurisdictional rapid habitat assessment data, and (2) trends in channel dimensions and hydraulics at USGS streamgages, referred to as specific gage analysis (Cashman and others, 2021), a methodological update to specific-stage analysis originally defined by Gilbert (1917). These two methods have different but complementary hydromorphic foci and provide different spatial coverages across the Chesapeake Bay watershed. Data compilation for specific-gage data and quality assurance and quality control (QA/QC) for rapid habitat data were completed in R (ver. 4.2.1; R Core Team, 2022b).

2.6.2. Analysis

Generalized additive models (GAMs) were used to evaluate trends across the Chesapeake Bay watershed in (1) multi-jurisdictional rapid habitat assessment data and (2) channel dimensions and hydraulics at USGS gages via specific gage analysis. Status was only evaluated for the multi-jurisdictional rapid habitat data. All analyses for both the specific-gage data and rapid-habitat data were completed in R (ver. 4.2.1; R Core Team, 2022b).

2.6.3. Results and Discussion

Results of GAMs for (1) multi-jurisdictional rapid habitat assessment data and (2) channel dimensions and hydraulics at USGS gages via specific gage analysis are conveyed and discussed. Results of the rapid habitat assessment analysis found evidence suggesting degrading hydromorphic conditions at the majority of sites with significant trends. However, the limited number of trend sites and spatial clustering of available sites prevented conclusions about hydromorphic conditions on a watershed-wide scale. Specific gage analysis revealed patterns consistent with channel evolution following disturbance, including decreasing bed elevation and increasing channel area trends indicating channel incision and expansion over a 75-year interval followed by occasional fill and widening, indicated by increasing bed elevation and channel area trends in the 10-year interval.

2.6.3.1. Rapid Habitat Assessment

On the scale of 0 (“poor” habitat quality) to 20 (“optimal” habitat quality), the status of rapid habitat assessment metrics averaged from 12.0 to 16.4. The lowest average scores were in Embeddedness (12.0), and the greatest average scores were in Channel Alteration (16.4). Status varied across the watershed. Low scores for some metrics were notably clustered in the Delmarva Peninsula and in the urbanized Baltimore–Washington, D.C., metropolitan area, and metrics with greater scores were generally in western parts of the watershed (fig. 20). The sites used for the analysis were also distributed unevenly across the watershed and tended to concentrate within Maryland and the southwestern portion of the watershed (Virginia and West Virginia).

Most sites that showed a change in habitat conditions were degrading, with relatively few sites showing positive, improving conditions. Of the 955 unique metric and site combinations tested for trends, 86 (9 percent) had “strong” trends, 96 (10 percent) had “possible” trends, 764 (80 percent) had “unlikely” trends, and 9 sites had no measured change across their record. Of those with “strong” trends, 54 (77 percent) were negative and indicated degrading condition, 13 (19 percent) were positive and indicated improving condition, and 3 results had nonmonotonic trends that indicated no net change in condition across the time interval. Although some metrics showed slightly different proportions of negative or positive trends, most metrics showed slight to moderate negative trends of degrading condition through time at most sites, and few sites showed positive, improving conditions (fig. 21). On average, linearized slopes for ”strong” trends changed by an absolute magnitude of change of about 0.55 points per year in either positive or negative directions, and the mean change across all metrics, –0.37 points per year, indicating overall degrading conditions. Across the trend interval of 2008–17, this would result in a change of nearly 5 points of change, and mean degradation of nearly 3.7 points, out of the 20-point scale, across our period of analysis. Areas with significant changes were scattered throughout the watershed, although with notable significant trends around the Baltimore–Washington, D.C., metropolitan area for many metrics (fig. 22).

Trend slopes were compared against metric scores from the beginning of the record to evaluate whether there were any management-relevant patterns in trend slopes based on a site’s initial condition. Epifaunal Substrate/Available Cover exhibited a significant relationship between trend slope and initial metric score, indicating that negative trends in epifaunal substrate were more likely to be observed at higher-quality sites (fig. 23). However, this trend analysis excluded some earlier data available at each site, which may provide additional perspective on how these sites may have changed through time. To achieve Chesapeake Bay Program targets, conserving sites currently in good condition and targeting lower quality sites for management interventions and (or) restoration may be necessary.

Overall, the field-based rapid habitat assessment data indicated that metrics of physical habitat condition that were significantly changing through time were degrading in condition, and only a few sites showed improvement. However, site locations included in this analysis were not evenly distributed throughout the watershed. Sites were instead clustered within a few jurisdictions and overrepresented many of the watershed’s metropolitan areas, which contain more urban stressors of physical habitat. Future analyses, which can incorporate additional data at more locations, may provide additional spatial and temporal perspective for these trends across the watershed. For example, if agencies and jurisdictions continue to collect rapid habitat assessments, additional sites should have enough data to meet the criteria for a trend analysis. Although representative of the conditions in these monitored rivers and streams, these locations might overlook other key drivers of change throughout the watershed, including areas of forest conservation, forest harvest, or mixed farming. This method also relies on visual-based scoring and subjective judgement, so uncertainty can be introduced by differing field personnel, differing standards among jurisdictions conducting the surveys, or guidelines to differentiate between similar scores in a similar quality category. Refer to section 2.8 “Indicator Results Synthesis” for more information.

2.6.3.2. Specific Gage Analyses

Across the 342 sites analyzed and the four trend intervals, all sites had at least one strong field measurement model that was sufficient for trend analysis. Overall, between 36 and 51 percent of trends analyzed in each time interval were significant at p-values less than or equal to 0.05. The greatest proportion of significant trends were within the 75-year trend interval. Among metrics, channel velocity and channel area had the least number of significant trends (41 percent), and bed elevation had the most (47 percent). Overall, trends were more likely to be significantly decreasing across the time interval for bed elevation and channel velocity, and the proportion of decreasing trends was greatest at 50- and 75-year intervals. In contrast, trends were more likely to be increasing for channel area, and the proportion of increasing trends was greatest at 50- and 75-year intervals (table 14). Overall, these patterns demonstrate the sequence of channel evolution following disturbance, like that from urban development, and demonstrate multiple channel evolution stages, including initial channel incision (negative bed elevation trends), expansion (increasing channel area), followed by occasional fill and widening (positive bed elevation trends and increasing channel area), and approaching new quasi-equilibrium conditions (indicating lesser change in metrics in recent years compared to previous change). Channel fill and constriction (decreasing channel area and increasing bed elevation) are also channel responses expected from increased upstream sediment supply. This also likely indicates decadal-scale channel adjustment to changing land use, rapid urbanization, and other development throughout regions of the Chesapeake Bay watershed because altered patterns of flow and sediment delivery have resulted in long-term changes and adjustments (Wolman, 1967). Widespread channel changes captured in 50- and 75-year trends (fig. 24) are likely driven by suburban expansion and other widespread change. More recent trend patterns, as evidenced in the 10-year trend interval (fig. 25), represent a mixed combination of recovery from former disturbance, more localized recent disturbances or development, and other site-specific factors, such as dam removals (Cashman and others, 2021).

Table 14.    

Number and direction of significant trend results in the Chesapeake Bay watershed for hydromorphology-specific gage analyses broken down by field metric and trend interval.

[Data are from Boyle and others (2025). Qualifying streamgage counts are in table 13. Trend interval years are 10 (water years [WY] 2012–22), 25 (WY 1997–2022), 50 (WY 1972–2022), and 75 (WY 1947–2022)]

Table 14.    Number and direction of significant trend results in the Chesapeake Bay watershed for hydromorphology-specific gage analyses broken down by field metric and trend interval.

Field metric	Trend interval, in years	Number of significant trends	Trend direction
			Percent of total decreasing	Percent of total increasing
Bed elevation	10	142	54.9	45.1
	25	101	57.4	42.6
	50	56	66.1	33.9
	75	30	73.3	26.7
Channel area	10	136	47.8	52.2
	25	76	38.2	61.8
	50	52	23.1	76.9
	75	28	17.9	82.1
Channel velocity	10	129	51.9	48.1
	25	79	58.2	41.8
	50	53	81.1	18.9
	75	33	84.8	15.2

This study empirically demonstrates the decadal-scale response of channel evolution and geomorphic change across many rivers throughout a region. In addition, this analysis highlights additional factors associated with changes in hydraulic character, such as flow-normalized velocity (refer to appendix 3), which co-occurs with channel evolution and geomorphic disturbance. As channel area increased in many rivers across the watershed, especially in response to urbanization, the flow-normalized channel velocity decreased, a habitat condition which may be relevant to many species, particularly lithophilic and other flow-sensitive species.

The trend results from specific gage analysis provide a comprehensive, watershed-wide approach for hydromorphic trends and multiple results covering jurisdictions throughout the Chesapeake Bay watershed. These trends highlight large-scale response to changes in watershed conditions over time but also locally driven patterns of geomorphic change demonstrating the timing and type of local watershed factors. As a result, individual streamgage results should be interpreted within the context of a local watershed’s history of disturbance and other changes. In addition, many USGS streamgages are operated to address specific management questions, such as a targeted restoration. These gages may provide additional local interpretive context but may not be generally representative of all streams in a region. Lastly, USGS streamgages are generally known to be biased toward larger rivers compared to the overall drainage network (refer to section 2.8 “Results Synthesis”); therefore, these results may not generally capture the range of channel changes occurring within headwaters. However, the results of the specific gage analysis capture unique information about channel change throughout the watershed and can help interpret the magnitude of geomorphic change and the timeframe for geomorphic recovery after disturbance (Cashman and others, 2021).

2.7.1. Data Compilation

Data were gathered from two Chesapeake Bay watershed-wide compilation efforts: Smith and others’ (2017) generation of the Chessie BIBI, which compiled benthic macroinvertebrate (hereafter termed macroinvertebrate) data; and Maloney and Krause’s (2021) development of a fish habitat assessment, which compiled fish sampling data. The two resulting datasets allow numerous approaches for which site-specific measures of benthic macroinvertebrates and fishes may be characterized, including multi-metric composite indices and summary metrics that describe assemblage composition (for example, life histories and tolerance classifications), taxonomic composition, and presence or abundance of specific species.

Long-term biological monitoring sites are generally sampled only once per year or less. To identify trend sites, datasets were separated by season [winter (December–February), spring (March–May), summer (June–August), and fall (September–November)], because season can have a strong effect on macroinvertebrate and fish community structure (Gelwick, 1990; Linke and others, 1999). To incorporate the largest number of sites possible within a trend interval long enough for trend analysis, the trend interval was set from 2008 to 2017. Sites qualified for trend analysis if samples were collected in the same season for at least 7 total years within the interval. The total needed to include at least one sample each in 2008 and 2017. If a trend site had multiple samples within the same season and year, one sample was randomly selected for analysis and the others were discarded. Following these steps, there were 92 spring, 4 summer, and 2 fall trend sites for macroinvertebrates and 36 summer and 8 fall sites for fish trend analyses (table 16).

2.7.2. Analysis

All analyses were completed in R (ver. 4.2.1; R Core Team, 2022b). Status of biological sites was defined as the average value of the 3 most recent years of data for each metric at each site. Both multi-metric indices, the macroinvertebrate Chessie BIBI and the fish MMI, have built-in narrative categories based on score ranges. The Chessie BIBI scores sites amongst five categories ranging from “very poor” to “excellent.” The fish MMI categorizes sites as “poor” or “good,” and sites in between these two categories are listed as “fair.” After averaging the index scores from the most recent 3 years at each site, the resulting values were placed into their corresponding narrative categories. The remaining three metrics (assemblage sensitivity, functional feeding strategy, and habitat preference) are presented as the percentage of the total assemblage composition, and the values from the last 3 years at each site were averaged to provide status.

For each site, trends were estimated separately for each macroinvertebrate and fish metric using GAMs. In total, 572 individual GAMs (396 for macroinvertebrates and 176 for fish) were fit using the R package “mgcv” (ver. 1.9-1; Wood, 2017). Because samples sizes were small (7–10 years per site with one measurement per year), the flexibility of the smooth function was limited for all GAMs by setting the upper limit on the degrees of freedom (k) associated with a smooth to be an integer no more than half of the total data points per trend. This allowed greater flexibility in model fit for sites with larger samples sizes but limited trend detection to highly linear trends for sites with small samples. Figure 26 displays the difference in function flexibility between a site with 7 data points and a site with 10 data points. Following methods from Murphy and others (2019), trends from each GAM were compared to a null model (GAM fit without year covariate) and were considered “strong” and “extremely likely” if the AIC value was at least 7 points lower than the null model and the p-value was less than or equal to 0.05. Results from this analysis can be found in Boyle and others (2025).

2.7.3. Results and Discussion

Patterns in status metrics varied greatly across the watershed but generally showed that around half of macroinvertebrate sites and one-third of fish sites were in poor condition and had few sensitive taxa and few taxa preferring swift water habitat. Analysis of status revealed that 56 of the 99 macroinvertebrate sites (57 percent) were in the “poor” or “very poor” Chessie BIBI condition categories, and 28 (28 percent) were in the “good” or “excellent” condition categories. At 21 sites (21 percent), less than 50 percent of the assemblage consisted of sensitive EPT-H taxa and clinger taxa (table 17; fig. 27). Sites in poor condition categories were clustered near the greater Baltimore–Washington, D.C., metropolitan area. Among the additional macroinvertebrate metrics, there were many more EPT-H and clinger taxa and fewer filterer taxa within assemblages in the Blue Ridge range of Virginia, compared to assemblages clustered near the Baltimore–Washington, D.C., metropolitan area (fig. 27).

Of the fish biological community sites analyzed, 14 (32 percent) were classified as “poor,” 9 (20 percent) were classified as “good,” and 21 (47 percent) had intermediate scores classified as “fair” (table 17; fig. 28). Percentages of nontolerant fish varied greatly among site locations. Percentages of rheophilic individuals followed geographic patterns of slope within the watershed. Assemblages generally consisted of less than 50-percent rheophilic individuals in sites near the Bay where low-slope gradients create slow water habitat, whereas assemblages in streams with swifter flows in the western portion of the watershed consisted of between 85 and 100 percent rheophilic individuals. Percentages of benthic invertivores varied greatly among sites with no consistent spatial pattern (fig. 28).

Many sites exhibited no strong directional trend in their biological metrics, indicating that biological composition quality is not increasing or decreasing at a detectable rate. Of the 572 trends in macroinvertebrate and fish metrics tested, 60 (11 percent) had strong increasing or decreasing trends (table 17). For both fish and macroinvertebrates, the assemblage sensitivity metric had the most strong trends of any metric, and most trends indicated declines in the percentage of sensitive organisms. Most sites had strong trends in only one metric (figs. 29, 30). Interestingly, strong trends in sensitive taxa did not often co-occur with strong trends in biotic index scores, highlighting the importance of analyzing individual metrics in addition to multimeric indices to capture patterns in specific groups of organisms that may be masked by an overall assemblage condition score. The lack of strong trends could be due to the small sample size and stringent analysis parameters rather than truly static biological condition through time. Though limiting the flexibility of the GAMs allowed for very high confidence in the reported strong trends, it likely prevented detection of trends with small sample sizes (7–8 years) and more variable or nonlinear patterns.

The current biological sites suitable for trend analysis are spatially biased; most sites are in Maryland, and some are also in and around the Blue Ridge range of Virginia. Sites with the most detectable trends were clustered around the highly urban Baltimore-Washington, D.C., metropolitan area or the protected and forested Blue Ridge Mountain (figs. 29, 30). There is no representation of the northern portion of the Chesapeake Bay watershed. Though the current dataset does not provide adequate coverage to interpret co-occurring trends or trend patterns for the entire Chesapeake Bay watershed, incorporating more data from a future data retrieval may help resolve both sample size and spatial distribution issues. Many sites have 5–6 years of data within the 2008–17 trend interval and are therefore close to meeting trend analysis criteria. If more recent data are identified for these sites, the overall number of sites suitable for trend analyses may greatly increase. Such recent data may also increase the sample size of current trend sites and the overall spread of trend sites across the watershed.

2.8.1. Trend Site Co-Occurrence

To assess the feasibility of a multi-indicator trend synthesis, trend sites of all seven indicators were overlayed to determine locations where multiple indicators co-occur. Sites representing multiple indicators that were clustered within a 200-meter radius of each other were identified as a potential co-occurrence location. Each cluster was visually inspected in ArcGIS Pro (ver. 3.3.2; Esri Inc., 2024) to ensure the potential cluster did not contain points falling on different streams. Most trend sites in this report were not co-located (table 18). There were no locations that had data for all seven indicators. One site, the Susquehanna River at Danville, Pennsylvania (USGS 01540500), had data for six of the seven indicators but did not have biological data (fig. 31). Indicators that regularly co-occurred were split into two general groups: (1) nutrients and sediment, salinity, flow, and stream temperature data, which were often collected at or near streamgages; and (2) biology and habitat data, which were often collected in the same sampling event by a stream sampling crew at wadeable locations (table 19). These two groupings highlight the mismatch between hydrology and water-quality data and biological community sample locations. Many gaged sites are not wadeable and difficult to accurately sample for biology. Similarly, many biological samples are collected at small streams, which are not well represented by the streamgaging network.

The limited instances of trend co-occurrence locations make integrating trend patterns across indicators difficult. Even when long-term data for multiple indicators co-occur spatially, differing trend intervals among indicators may prevent trend comparison unless enough data exists to truncate datasets to uniform start and end dates. Integrating trend results from multiple indicators at specific locations will require new and creative methods.

2.8.2. Status Snapshot

Although synthesizing trends results across indicators remained a challenge, status conditions were synthesized across indicators, and potential drivers were explored by assessing status conditions during a period shared among all indicators. The most recent status period shared by all indicators (calendar years 2015 through 2017) was selected as the common period of interest to assess if land cover is a driver of status condition. The salinity, hydromorphology, and biological assemblage indicators had status values that already corresponded with this time interval. For the status snapshot, status values for salinity departure from background SC, all nine hydromorphological metrics, and the multi-metric indices and assemblage sensitivity metrics for fish and macroinvertebrate biological assemblages were utilized. For the remaining indicators (stream temperature, nutrients and sediment, and streamflow), additional status metrics were calculated from 2015–17 data for each site. Temperature was represented by mean temperature, nutrients and sediment were represented by mean nutrient and suspended sediment concentrations, and streamflow was represented as the mean flow metric values for 2015–17.

The spatial distribution of the highest and lowest quality sites for each indicator was assessed using data from the shared status interval for all indicators. High-quality sites were sites with the status metric values most indicative of good stream condition for each indicator, whereas low-quality sites represent the poorest status values for each indicator. Metric values were divided into quartiles and coded as “degraded condition” or “good condition” if the status value fell within the top or bottom 25 percent of values. Which percentile was coded as degraded or good varied among the metrics. For example, 75th-percentile temperature values were categorized as degraded condition, whereas 75th-percentile macroinvertebrate biological index scores were categorized as good condition. To categorize flow metric values as good condition or degraded condition, the percentage difference between mean annual flow metric values for 2015–17 and the mean annual metric values for 1985–2022 (the full interval of flow data) was calculated and divided into quartiles.

After coding, sites for each metric were then assigned to the overall high-, low-, or intermediate-quality categories. Temperature sites were assigned the quality category that corresponded to their single temperature status metric (mean temperature). For example, sites with degraded temperature values were assigned to the low-quality temperature site category. For the nutrients and suspended sediment and biology indicators, a low-quality site was one associated with any degraded metric(s), and a high-quality site was one associated with any good metric(s) and no degraded condition metric(s). The streamflow and hydromorphology indicators had many metrics with sometimes conflicting scores. For these indicators, a site was categorized as low-quality if it was associated with any degraded metric(s) and no good metric(s), and vice versa for high-quality categorization. Sites with both good and degraded metrics were assigned to the intermediate category because of uncertainty. For salinity, a different definition was used; sites were considered low-quality if the status metric was greater than 3 times the predicted background SC and high-quality if the status value was at or below background SC. Status values and quality categories for each indicator site can be found in Boyle and others (2025).

The high-quality and low-quality sites were mapped individually to compare the overall spatial distribution of high- and low-quality sites. High-quality sites for all indicators were spread mostly evenly across the watershed, but some distinct clusters of high-quality biological assemblage, salinity, nutrients and sediment, and hydromorphology sites were apparent along the Blue Ridge range in Virginia and along the northern West Virginia–Maryland border (fig. 32). In contrast, the low-quality sites for multiple indicators were clearly clustered around the urban centers of Washington, D.C., and Baltimore, as well as within the Delmarva Peninsula (fig. 32). Maps showing the spatial distribution of high- and low-quality sites for each indicator individually can be found in appendix 4.

Potential drivers of status conditions were further investigated through a simple linear regression, which was used to determine if different major land-cover categories explained status values for each indicator. To characterize land cover for each status site, sites were associated with flowlines from the National Hydrography Dataset Plus High Resolution (NHDPlus HR) network and were subsequently joined to a dataset that summarized the total watershed area and National Land Cover Database (NLCD) land cover for each flowline (U.S. Geological Survey, 2022; Gressler and others, 2023;). The 2016 NLCD dataset was used because 2016 marked the middle of the 3-year common status interval. Multiple land-cover percentages were summed to create the following predictor variables: (1) the percentage of developed land cover (open space, low intensity, medium intensity, and high intensity), (2) the percentage of agriculture land cover (pasture or hay, and cultivated crops), and (3) the percentage of forest land cover (deciduous, evergreen, and mixed). Individual linear models were run for each combination of indicator metric and predictor variable and categorized relationships as significant at p-values less than or equal to 0.05. Linear model results are shown in table 20 and appendix 4 and can be found in Boyle and others (2025).

Many metrics representing the seven indicators had significant relationships with the land-cover predictor variables. Of the 24 metrics, 14 (58 percent) had significant relationships with developed land cover. High percentages of developed land cover were significantly associated with higher total phosphorus, stream temperature, and stream salinity; less extreme low-flow magnitudes; greater frequency of low and high-flow events; and shorter duration of high-flow events than low percentages of developed land cover. High percentages of developed land cover are also significantly associated with low values of habitat metrics, including bank stability, embeddedness, and epifaunal substrate, and low values of metrics characterizing macroinvertebrate sensitive taxa and assemblage condition. Fewer metrics (11 of the 24 metrics [46 percent]) had significant relationships with agriculture land cover. High percentages of agriculture land cover was associated with elevated total nitrogen, total phosphorus and, salinity concentrations; less extreme low-flow magnitudes; higher habitat velocity-depth combination scores; decreased high-flow duration, frequency, and magnitude; lower habitat streambed sedimentation scores; and fewer sensitive fishes (nontolerant individuals) compared to low percentages of agriculture land cover. Conversely, indicator responses to forest land cover were a near-perfect inverse of the responses to developed land cover. Forest land cover was significantly associated with lower nutrient and SC concentrations, lower stream temperatures, less frequent high-flow events, lower low-flow magnitude, longer high-flow duration, larger high-flow magnitude, and higher habitat scores and metrics for sensitive macroinvertebrates.

Together, the patterns between indicator status conditions and land-cover predictors demonstrate the negative impact of anthropogenic activity on stream condition. Urban development increases the amount of impervious surface and necessitates deicer application, leading to high salt loading into streams and groundwater systems (Kaushal and others, 2005). Increased sedimentation as a result of urban development also effects habitat quality. The removal of riparian buffers and alterations to natural flow regime that are prevalent in urban areas can cause additional degradation of instream habitat and increases in streamflow flashiness and stream temperature (LeBlanc and others, 1997; Paul and Meyer, 2001; Nelson and Palmer, 2007; Rosburg and others, 2017; Anim and others, 2018). High percentages of agriculture land cover were similarly associated to poorer water quality than watersheds with low percentages of agriculture land cover as agricultural practices can introduce high amounts of nutrients and increase SC in streams (Kaushal and others, 2018). Water withdrawals for agricultural practices can also alter flow regimes, creating more extreme low-flow events (Eheart and Tornil, 1999; Almahawis and others, 2024). Poorer water quality, flow, and habitat conditions associated with these anthropogenic activities were especially represented in the macroinvertebrate metrics. Overall macroinvertebrate assemblage condition scores (Chessie BIBI), and the percentage of sensitive taxa (EPT-H) were significantly lower at sites with greater urban development and the percentage of sensitive fish (percentage of nontolerant individuals) was significantly lower at sites with high percentages of agriculture land cover.

The association of indicators to the percentage of forest land cover was a near-perfect inverse relationship to that of developed and agriculture land covers. This is unsurprising, because the percentage of forest land cover is often negatively correlated with developed and agriculture land cover, but these results further highlight the importance of forested watershed and riparian buffers for maintaining good stream condition. Forested watersheds and riparian buffers reduce runoff of nutrients and shade streams, limiting stream warming (Lowrance and others, 1984; Bowler and others, 2012). Forested areas often have higher-quality habitat for biological assemblages by limiting streambank erosion and allowing for natural flow regime processes as shown by higher amounts of forested land covers being associated with less extreme low-flow magnitudes and less frequent high-flow events (Wang and others, 1997; Schoonover and others, 2006).

2.8.3. Land Cover and Spatial Representation of Indicator Site Networks

Multiple factors influence the location, frequency, and duration of the stream monitoring efforts highlighted in this study. Some sites may have been selected randomly, others may have been targeted for long-term sampling to monitor a pre- and (or) post-disturbance event (such as dam removal, bridge construction, or restoration project), and yet others may have been selected as longstanding reference sites in least-disturbed areas. To begin understanding how the indicator sites in this report represent the Chesapeake Bay watershed as a whole, cumulative-distribution plots were created for watershed area, urban, agriculture, and forest land covers for all indicator site networks (fig. 33). These distributions were compared visually to cumulative-distribution curves of all NHDPlus HR stream segments in the Chesapeake Bay watershed.

Small watershed settings were poorly represented by all indicator site networks except biological aquatic community sites. Biological aquatic community sites most closely followed the distribution of Chesapeake Bay watersheds because biological data are most often collected in the small, wadeable streams that make up most of the Chesapeake Bay watershed’s stream network. Toxic contaminants, nutrients and sediment, and salinity site networks had the highest proportions of large watersheds; 75 percent of sites in each network represented watersheds larger than 2,093 square miles (mi²), 928 mi², and 544 mi², respectively. In comparison, only 22 percent of Chesapeake Bay watershed NHDPlus HR segments represent watersheds greater than 40 mi². Areas with intensive agriculture (greater than 50 percent agriculture land cover) were under-represented by all indicator networks compared to agricultural intensity across all streams in the Chesapeake Bay watersheds. Distribution of urban settings across the salinity network was similar to the full watershed, but all other indicators represented more urban land cover within their watersheds than is represented across the full Chesapeake Bay watershed. All indicator sites had somewhat similar representation of forested land cover compared to the full Chesapeake Bay watershed, though the salinity network had more sites with moderately forested watersheds (50-percent forest land cover) compared to other indicators.

In addition to differences in land-cover representation compared to the full Chesapeake Bay watershed, indicator site networks have clustered and disparate spatial coverage across the watershed. Almost all indicators had low to no representation of sites in New York, except for streamflow, and all indicators had more sites in the southern half of the watershed compared to the northern half. Many indicator networks, especially temperature, stream habitat, and biology, had limited data for trend analysis, and criteria for trend inclusion were chosen to maximize the spatial coverage across the watershed and to provide the adequate record lengths for trend analysis. Many biological and stream habitat sites in the northern half of the watershed are close to meeting trend criteria, and temperature monitoring at gaged sites is becoming more prevalent. Thus, the spatial coverage and length of the trend intervals for all indicators will likely increase over time if data collection continues.