Introduction

Estuaries are currently facing significant challenges due to the influx of anthropogenic nutrients resulting from various human activities, including the growth of coastal populations, agricultural practices involving fertilizer and manure applications, and the introduction of nutrients through atmospheric deposition (Breitburg et al. 2018). Excess amounts of nutrients such as nitrogen and phosphorus can enhance algae growth and biomass, subsequently leading to eutrophication (i.e., organic matter production) and bottom water oxygen depletion (Paerl et al. 1998; Pinckney et al. 2001; Breitburg et al. 2018). A variety of long-term patterns in oxygen depletion have been observed across estuaries. Some systems exhibit persistent eutrophication and hypoxia due to consistently high nitrogen loading (e.g., Gulf of Mexico, Rabalais et al. 2001), while other systems have seen increasing nutrient loading and hypoxia extent with ongoing eutrophication (e.g., Bohai Sea, Xin et al. 2019). Alternatively, some systems such as the Chesapeake Bay and Long Island Sound have seen modest reductions in aspects of hypoxia after implementing nutrient reduction strategies (Murphy et al. 2011; Whitney and Vlahos 2021). In general, attempts to reduce hypoxia through nutrient management have been complicated by persistent nonpoint sources, increasing land use change, and high year-to-year variability in physical forcing (Rabotyagov et al. 2014).

Climate change can also alter hypoxia in estuaries through alterations in precipitation patterns, sea level rise, and warming (Diaz et al. 2019). Elevated precipitation can drive larger inputs of freshwater, which has the potential to increase vertical stratification and deliver more nutrients that can lead to an increased hypoxia (Testa et al. 2022). In some climates, reduced precipitation with climate change can reduce flushing rates and enhance hypoxia (Cottingham et al. 2018; Warwick et al. 2018). Bottom water hypoxia can also be intensified through sea level rise if it increases vertical stratification (Ni et al. 2019). Warming is expected to increase hypoxia because higher temperatures would decrease the solubility of oxygen in the water column (Weiss 1970) and potentially accelerate metabolic rates, such as nutrient recycling, phytoplankton productivity, and respiration to increase oxygen consumption (Najjar et al. 2009). Most projections of future climate in large estuaries suggest that warming will exacerbate hypoxia (Wang et al. 2017; Laurent et al. 2018; Irby et al. 2018; Ni et al. 2020; Meier et al. 2021), while other studies have shown that contemporary warming has already overcome expected declines associated with nutrient load reduction (Ni et al. 2020; Frankel et al. 2022).

Hypoxia can alter benthic community structure through direct and indirect pathways within the food web (Hale et al. 2016; Woodland et al. 2022). In particular, hypoxia can contribute to the mortality of fisheries and prey species, increase the susceptibility of organisms to disease, and reduce growth and production rates (Diaz et al. 2019). Dissolved oxygen response thresholds for benthic fauna typically begin around 2 mg L−1 with larger, epifaunal taxa such as many echinoderms and crustaceans generally more susceptible to the effects of lower oxygen levels than infauna such as mollusks (Vaquer-Sunyer and Duarte 2008; Levin et al. 2009). Mobile organisms will vacate hypoxic areas (Brady and Targett 2013), while organisms that remain in hypoxic waters are vulnerable to physiological stress or mortality, as well as behavioral changes that can lead to a temporary increase in predation on these organisms (Pihl et al. 1992; Long and Seitz 2008). Longer and more extensive hypoxic events increase production from microbial communities, reducing the energy transferred from primary to secondary trophic levels (Sturdivant et al. 2014). Understanding the development and drivers of hypoxia at multiple spatial scales over a long period of time will allow for the development of novel solutions to address linkages between hypoxia and food webs.

Hypoxia responses to nutrient loading and climate change have been extensively studied in large systems like Chesapeake Bay (Hagy et al. 2004; Testa and Kemp 2012; Li et al. 2016; Irby et al. 2018; Ni et al. 2020) and the Baltic Sea (Conley et al.2007, 2009; Meier et al. 2018, 2021; Savchuk 2018), while fewer studies have investigated smaller, shallower hypoxic zones. Thus, an improved understanding of the interactions of eutrophication and warming on hypoxia expansion (or reduction) in the wide variety of different estuarine environments is needed. The purpose of this paper is to quantify the long-term controls on hypoxia in a small, eutrophic estuarine basin to understand key biogeochemical and physical drivers and associate those long-term changes with benthic invertebrate biomass. We evaluated the Patuxent River estuary, which is a well-studied tributary of the Chesapeake Bay with oxygen depletion at multiple scales, a history of nitrogen loading reductions, and persistently degraded water quality (Kemp et al. 2005). We seek to answer questions that have not been previously addressed in detail for this estuary, including (1) How have the patterns of dissolved oxygen changed over time and space?; (2) How has hypoxic volume changed over time and what factors control it?; (3) How do the onset, duration, and breakup of hypoxia vary over time and space?; and (4) Is there a relationship between hypoxia and benthic biomass? We examined long-term, publicly available datasets involving water column biogeochemical variables, riverine discharge, and physical factors combined with estimates of watershed nutrient loading. Our findings advance knowledge around the interaction of eutrophication and climate change in controlling oxygen depletion in coastal waters and the ecological linkages between hypoxia and benthic communities.

Methods

Description of Study Stations and System

The Patuxent River is one of the tributaries of the Chesapeake Bay that exhibits a two-layered circulation pattern with a landward-flowing bottom water layer and a seaward-flowing surface layer (Testa et al. 2008). The river has a maximum depth of 53 m and is a major social, environmental, and economic resource for the region (Voinov et al. 1999; Costanza et al. 2002; Bartell 2003; Breitburg et al. 2003b; Fisher et al. 2006). Following upgrades of sewage treatment plants, the tidal fresh regions of the Patuxent River estuary have shown improvements in water quality, resulting in more oxygen availability and declines in bottom-layer oxygen consumption; however, the lower estuary water quality has been degraded since the mid-1990s due to an increase in riverine discharge and dissolved inorganic nitrogen from the Chesapeake Bay (Testa et al. 2008). For this study, the estuary was examined across three distinct zones: (1) the tidal fresh (TF) portion that is shallow (< 5 m), low salinity, and typically vertically well-mixed; (2) the riverine-estuarine transition (RET) portion that is turbid and mesohaline; and (3) the lower estuarine (LE) portion that includes a deep (> 10 m) central channel and exchange is at the seaward end with the Chesapeake Bay (Fig. 1). We examined data from eight long-term monitoring stations, including stations TF1.5 and TF1.6 that can be categorized as freshwater (salinity = 0–2 PSU); stations TF1.7 and RET1.1 that can be categorized as lower mesohaline (salinity = 5–10 PSU); and stations LE1.1, LE1.2, LE1.3, and LE1.4 that can be categorized as higher mesohaline (salinity = 10–16 PSU). We will present time series data from 1985 to 2021 and explore spatial patterns at representative stations TF1.6, RET1.1, and LE1.1.

Fig. 1
figure 1

Map of the Patuxent River estuary, including location of water quality (squares) and benthic invertebrate (blue circles) monitoring stations (TF, tidal fresh; RET, riverine-estuarine transition; LE, lower estuary)

Full size image

Data Sources

Vertical profiles of dissolved oxygen (DO, mg L−1), salinity, and water temperature (°C) from 01/01/1985 to 12/31/2021 were obtained from the Chesapeake Bay Data Hub (www.chesapeakebay.net) for the following stations: TF1.5, TF1.6, TF1.7, RET1.1, LE1.1, LE1.2, LE1.3, LE1.4. The data were collected by the Maryland Department of Natural Resources 1–2 times per month as part of the Chesapeake Bay Program’s Water Quality Monitoring Program. From 1985 to 2009, measurements were taken twice a month with the exception of November–February, while from 2010 to 2021, measurements were taken monthly. We analyzed surface and bottom DO, temperature, and salinity measurements, where surface measurements were made at 0.5 m from the surface and bottom concentrations only included the lowest recorded depth on a given sampling day. Measurements were taken with Yellow Springs Instrument (YSI) data sondes and HydroLab multiparameter sondes.

Estimates of watershed nutrient (nitrogen, phosphorus) and freshwater input rates over 1985–2021 were obtained from a combination of watershed model estimates and statistical model estimates from a United States Geological Survey (USGS) monitoring station at Bowie, MD (01594440), using Weighted Regressions on Time, Discharge, and Season (WRTDS) analysis (Zhang et al. 2022; Mason and Soroka 2022). The watershed nutrient loading and flow rates outside of the gauged regions were estimated by the Chesapeake Bay Phase 6 model (Chesapeake Bay Program 2020), which generated predictions for all loads below the fall line (BFL), including nonpoint source loads below the fall line (nps.BFL), point source loads below the fall line (ps.BFL), and shoreline erosion (sho.BFL). Loading estimates were made for specific forms of nitrogen and phosphorus, including dissolved inorganic nitrogen (DIN), total nitrogen (TN), dissolved inorganic phosphorus (DIP), and total phosphorus (TP).

To quantify long-term changes, the average and standard deviation of DO, temperature, salinity, and chlorophyll-a (chl-a) were computed for surface and bottom water of each station for each year. We calculated annual averages and seasonal averages, where summer was defined as June–August, fall as September–November, winter as December–February, and spring as March–May. To isolate the biological vs. physical drivers of the DO time series, we calculated apparent oxygen utilization (AOU = DOsat – DOobs), where DOobs is the measured value and DOsat is the oxygen saturation concentration determined at in situ temperature and salinity following Weiss (1970).

Hypoxic Volume and Hypoxic Volume Days (HVD)

We interpolated DO profile data for all available dates in the dataset to a two-dimensional grid (1 m (vertical) by 1.8 km (horizontal)) using natural, linear, and cubic interpolation methods (Figs. S2.1–S2.3). The grid spans 56 km horizontally along the estuarine longitudinal axis and up to 34 m deep. The volumes of the water within the interpolated 2D grid were then computed with axial depth and volume data (Cronin and Pritchard 1975). Hypoxic volume was calculated for each sampling period by summing the volume of all interpolated cells with an oxygen concentration of less than 2 mg L−1. We then computed HVD for each year by linearly interpolating the estimated hypoxic volumes to daily values and then taking the sum of the daily estimated hypoxic volume time series for each year. To quantify the variance in HVD over time, z-scores were calculated across two time periods, 1985 to 2010 and after 2010.

We compared our estimates of HVD to two alternative interpolation techniques. First, we compared our estimates to those of Testa et al. (2008), which used a linear interpolation of the oxygen concentrations and the same 2D grid as we used in our estimates. These values were similar to our estimate (5.05 ± 3.03 km3) and were highly correlated (r2 = 0.96). We also compared our estimates with those obtained from an inverse distance weighting (IDW) approach used by the Chesapeake Bay Program, which computes 2D/3D interpolation for the Bay and its tributaries (Bahner 2006; Fig. S2.4). The method involves cells that are 1 km east–west by 1 km north–south and 1 m vertical from surface to bottom, and each cell midpoint is computed from the four nearest neighboring water quality values (Bahner 2006). These estimates were ~ 40% lower than our estimates (2.68 ± 1.87 km3), but highly correlated (r2 = 0.90). We determined that the differences between the two methods were correlated to the magnitude of the HVD (r2 = 0.81, p < 0.001), suggesting that our approach had a possible positive bias under conditions of larger volumes. We tested the impact of this bias on our HVD estimates by removing this bias from the HVD estimate (HVDcorr = HVD – 0.47*HVD), where HVD is our estimate and 0.47 is the slope of the relationship between our HVD estimate and the difference between our HVD and the IDW-estimated HVD (Figs. S2.4–S2.5). We determined that the upward trend in our HVD data was maintained following this correction (r2 = 0.11, p = 0.048), albeit with a slightly smaller slope (y = 0.06*HVD – 112). Thus, the temporal trend in our HVD estimates should not be biased by our choice of method. We also computed hypoxia duration for the entire estuary by summing the total number of days in the HVD time series when hypoxia occurred and also by computing the difference between the first-day hypoxia occurred in a given year and the last day (see Figs. S2.6–S2.7 for details).

Oxygen Depletion Rates and Hypoxia Onset Dates

Oxygen depletion rates were calculated for each station and each year from 1985 to 2021 during the months of March–July by applying a linear regression of the observed temporal decline in bottom water oxygen concentration (see Figs. S3.2–S3.9 for detailed data for each year and station). The slope was extracted from these models and used as the oxygen depletion rate (mg L−1 d−1). Hypoxia onset dates were found by first taking a cubic interpolation between the monthly data using the spline method for each station to create daily DO concentrations (Testa and Kemp 2014; Fig. S3.1). Hypoxia onset and breakup dates, defined as the day of the calendar year in which DO fell below 2 mg L−1 and equaled or was greater than 2 mg L−1, respectively, were generated for each station for each year using the interpolated values. To test for differences between stations, a one-way ANOVA followed by a Tukey’s HSD test was performed using R Statistical Software (v4.2.1; R Core Team 2022). A sensitivity analysis was run to ensure that the number of sampling times in a month (1–2) did not bias the results. This entailed pulling one sampling record from each month and rerunning the analyses and then plotting a linear regression to see the correlation between the full results vs. the results with one sampling day selected. The sensitivity analysis showed that the number of sampling times in a month did not alter the interannual variability or timing of onset (r2 = 0.76; p = 0.55). However, the timing of breakup was affected (r2 = 0.54, p = 0.02). Breakup dates tended to be earlier when there were two sampling events in a month.

Stratification Strength

To gain insights into the vertical structure of the water column for each year and its respective season, we calculated an estimate of stratification strength (ΔσT) for all stations and seasons by subtracting the density of bottom water from the density of surface water. Water density (σT) was calculated using the in situ temperature and salinity conditions, following Lewis (1980). If the stratification strength was greater than 1.5, the water column was categorized as “stratified,” indicating a significant density difference between surface and bottom waters. On the other hand, if the stratification strength was less than 1.5, the water column was classified as “well-mixed,” implying a relatively homogeneous density distribution throughout the water column (Levitus 1982).

Watershed Input Rates and Temporal Trends

The non-parametric Mann-Kendall tests were conducted to quantify the annual directionality, strength, and significance of temporal trends in TP (load), TN (load), and Q (discharge) to the estuary from both the USGS gauging station (hereafter RIM = river input monitoring) and BFL loads. We used Pearson correlation (rP = correlation coefficient) to compare variations in HVD, TP, TN, discharge, TNnps.BFL, TPnps.BFL, and Qnps.BFL. We also used linear regression to compute temporal trends in the time series of all water column variables including bottom water DO, temperature, salinity, AOU, stratification, HVD, oxygen depletion rate, hypoxia onset date, surface and bottom water chl-a, benthic biomass, and watershed inputs of discharge, TN, and TP.

Benthic Invertebrate Biomass

We also estimated benthic biomass from data collected at several long-term monitoring stations for benthic invertebrates. The monitoring of benthic species biomass in the Maryland Chesapeake Bay mainstem and tributaries has been carried out since July 1984 as part of the Chesapeake Bay Benthic Monitoring Program (2023). The dataset includes sampling from fixed locations in May and August–September, as well as random locations in August–September each year. In the Patuxent estuary, data were collected from three fixed stations (71, 74, and 77) in May (spring) from 1985 to 2008 and in August–September (summer) from 1985 to 2021 (Fig. 1). The fixed sites were sampled twice a year through 2008, in May and in late August or September. Subsequently, starting from 2009, the sites have been sampled once a year in late August or September. To ensure consistency, we considered both August and September data as part of the summer season. At each of these fixed sites, three replicate sediment samples for benthos were collected using the same gear employed since 1984. At these three stations, the collection method for ash-free dry weight in grams includes the use of the 220 cm2 surface area Wildco box corer benthic grab. The contents obtained were sieved through a 0.5-mm screen, thus targeting macrofauna, and preserved in the field in a 10% formaldehyde solution. We kept the biomass values in grams and averaged the data across all three replicates to obtain mean values per station. Detailed information about the collection methods and study design can be found in the Data Dictionary for Maryland and Virginia (https://www.baybenthos.versar.com/data.htm), ensuring transparency and accuracy in the data analysis.

We divided the benthic data into summer and spring sets. Given the highly skewed distribution of benthic biomass data, we applied loge-transformation, obtaining the natural log of grams dry weight biomass. Using the transformed data, we fitted linear regressions to HVD (predictor) and the summer and spring loge-transformed biomass data (response variables). To assess the relationships between HVD and benthos, we examined two scenarios: one incorporating the total biomass, averaged across three replicates and considering all species, and another focusing on three distinct taxonomic groups of benthic invertebrates. These taxonomic groups included all species belonging to Bivalvia (hereafter, “bivalves”), Malacostraca (“small crustaceans,” includes amphipods, isopoda, mysids, tanaids, and cumaceans), and a combined Annelida and Nemertea group (“worms,” includes polychaetes, oligochaetes, and nemerteans), with detailed information on species composition provided in Table S1. The regression analysis aimed to capture the relationships between HVD and each functional group’s biomass during both summer and spring.

Results

Temporal and Spatial Trends in DO, Temperature, and Salinity

Surface water oxygen concentrations did not increase or decrease over time at any station, and the highest DO concentrations occurred during winter (Fig. 2a–c). In the tidal fresh portion of the river, represented by station TF1.6, bottom water DO decreased in fall and summer over time, with no significant changes in winter and spring (Fig. 2d–f; Table 1a). In the riverine-estuarine transition portion, represented by RET1.1, bottom water DO also decreased in summer and fall, while no considerable changes occurred in winter and spring. In the lower estuarine portion, represented by LE1.1, bottom water DO decreased in winter and fall and did not change in spring and summer. Only the tidal fresh portion of the river shows significant decreases in annual mean bottom water DO concentrations (Table 1a). We found an increase in temperature in all seasons but spring (Table 1; Fig. S5.1) and no trends in salinity (Table 1; Fig. S5.2).

Fig. 2
figure 2

Mean DO concentration time series for each season from 1985 to 2021, spanning the tidal freshwater seaward to the mesohaline region. Summer was defined from June–August, fall from September–November, winter from December–February, and spring from March–May. Surface (a–c) and bottom concentrations (d–f) are included for each of three representative stations along the estuarine axis. Total depths for each station are 6 m (TF1.6), 11 m (RET1.1), and 12 m (LE1.1)

Full size image
Table 1 Trends in computed metrics of oxygen availability and related physical and biogeochemical variables computed over 1985–2021 in the tidal fresh (a), riverine-estuarine transition (b), and lower estuarine zones (c) in the Patuxent River estuary. Trends were derived from linear regression of the seasonal- or annual-scale metrics. Increases (p < 0.1) are denoted with an upward arrow, decreases with a downward arrow (p < 0.1), and no change (p > 0.1) with a straight line. In regard to water quality, green shading indicates a positive trend relative to oxygen levels, red shading indicates a negative trend relative to oxygen levels, and yellow indicates no trend. N/A indicates that a trend was not able to be computed for that time scale. The sample size for each dataset is indicated by the following symbols: Δ for n = 12, ♦ for n = 25, * for n = 37, & for n = 36, % for n = 35, and $ for n = 34. Bottom chlorophyll-a time series can be found in Figs. S5.4 to S5.7
Full size table

Hypoxic Volume Days

Computed HVD increased over the 1985–2021 period, at a rate of 0.11 km3-d y−1 (r2 = 0.13, p = 0.03; Fig. 3a). HVD variability also increased over time (Fig. S5.3) with a z-score of -0.13 before 2010 and 0.43 after 2010. The increasing variance after 2010 was associated with the two highest HVDs in 2011 and 2018 and the two lowest HVDs in 2012 and 2017, illustrating that consecutive years can be dramatically different.

Fig. 3
figure 3

Estimates of HVD over 1985–2021 (a) computed from daily interpolations of monthly estimates of hypoxic volume (b; see methods for details on interpolations). The black trend line in the top panel is the linear regression (n = 37) of the annual estimates of HVD against time

Full size image

Oxygen Depletion Rate and Hypoxia Duration

The duration of hypoxia in the Patuxent River estuary ranged from 0 days in 2012 to 194 days in 2020 (Fig. 4). When dividing the entire time span into three time periods, the average duration of hypoxia in the later time period (2010–2021) exceeded that of the early period (1985–1996) by 13.6 days. Given that no hypoxic days were recorded in 2012 and if the zero value in 2012 is disregarded, the difference in the average number of hypoxic days between the two periods increased to 24.7 days. Hypoxia duration also varied spatially in the Patuxent River estuary. Hypoxia tended to originate at station LE1.1 at the beginning of June and spread up-estuary toward RET1.1 and down-estuary toward LE1.2 by mid-June (Fig. 5a). LE1.1 also stayed hypoxic longer than the other stations and tended to remain hypoxic until mid-August. Stations RET1.1, LE1.2, LE1.3, and LE1.4 had breakup days occurring in mid-July, with hypoxia lasting for about a month. There was no trend in hypoxia onset or breakup date from 1985 to 2021 (Table 1, Fig. S3.1). Pearson correlations between onset day and winter- and spring-averaged bottom water chl-a, stratification strength, and river flow did not explain substantial amounts of variability in onset day across the estuary (Figs. S4.1–S4.12). The rate of water column oxygen depletion varied spatially in the Patuxent River estuary (Fig. 5b). The LE sites had the highest depletion rates (0.08 ± 0.01 mg L−1 d−1) over the full time series, and the TF stations had significantly lower oxygen depletion rates (p < 0.05) and did not undergo hypoxia (Fig. 5b). The oxygen depletion rate at RET1.1 was higher than the TF stations and lower than the LE stations (Fig. 5a). There was no trend in oxygen depletion rate from 1985 to 2021 at representative stations TF1.6 (p = 0.76), RET1.1 (p = 0.45), or LE1.1 (p = 0.43) (Table 1).

Fig. 4
figure 4

Patuxent River estuary annual hypoxia (DO < 2 mg L−1) duration in 1985–2021, estimated from counting the number of days that hypoxia was recorded in each year. The numbers displayed above the bars indicate the count of hypoxic days for each year, with no recorded hypoxia in 2012. The 37-year period is segmented into three intervals, each marked by its average and coefficient of variation: 1985–1996, 1997–2009, and 2010–2021

Full size image
Fig. 5
figure 5

Box plots of hypoxia onset/breakup date (a) and oxygen depletion rate (b) for eight stations in the Patuxent River estuary, spanning the tidal freshwater region (e.g., TF1.5) seaward to the mesohaline region (e.g., LE1.4). Left boxes are hypoxia onset dates and right boxes are hypoxia breakup dates for each station (a). Box plots are based on 1985–2021 estimates and include long-term median values (center line), 25th and 75th percentile (bottom and top of the boxes, respectively), the range between minimum and maximum values (vertical bars), and outliers (circles). The color scheme for water quality monitoring stations corresponds to that in Fig. 1. Total depths for each station are 10 m (TF1.5), 6 m (TF1.6), 3 m (TF1.7), 11 m (RET1.1), 12 m (LE1.1), 17 m (LE1.2), 23 m (LE1.3), and 15 m (LE1.4)

Full size image

Apparent Oxygen Utilization

The calculation of AOU revealed a consistent pattern across seasons in the surface waters of TF1.6, RET1.1, and LE1.1 (Fig. 6a–c), while bottom waters exhibited more seasonal variability in AOU (Fig. 6d–f). The bottom waters showed increasing AOU in fall across all regions of the river (p < 0.001) and did not show trends for other seasons (Fig. 6, Table 1).

Fig. 6
figure 6

Mean AOU time series for each season from 1985 to 2021, spanning the tidal freshwater seaward to the mesohaline region. Surface AOU (a–c) and bottom AOU (d–f) are included for each of three representative stations along the estuarine axis. Total depths for each station are 6 m (TF1.6), 11 m (RET1.1), and 12 m (LE1.1)

Full size image

Stratification Strength

There was minimal stratification in the tidal fresh region of the estuary (TF1.5, TF1.6, TF1.7), but substantial vertical stratification in the middle of the estuary (RET1.1 and LE1.1), giving way to modest stratification in seaward regions (LE1.2, LE1.3, LE1.4; Fig. 7). In the middle and lower part of the estuary (e.g., RET1.1, LE1.1), stratification strength increased over the record during fall and annually (Table 1, Fig. S5.3), but there were no trends in stratification strength in the well-mixed tidal fresh estuary.

Fig. 7
figure 7

Box plots of seasonal stratification strength at eight stations in the Patuxent estuary spanning the tidal freshwater region (e.g., TF1.5) seaward to the mesohaline region (e.g., LE1.4) in spring (a), summer (b), fall (c), and winter (d). Box plots are based on 1985–2021 estimates and include long-term median values (center line), 25th and 75th percentile (bottom and top of the boxes, respectively), the range between minimum and maximum values (vertical bars), and outliers (circles). Total depths for each station are 10 m (TF1.5), 6 m (TF1.6), 3 m (TF1.7), 11 m (RET1.1), 12 m (LE1.1), 17 m (LE1.2), 23 m (LE1.3), and 15 m (LE1.4). The horizontal red line is the index we used to identify stratification

Full size image

Watershed Input Rates and Temporal Trends

TP loads varied over the time series from 10 to 750 kg d−1, TN load varied from 1000 to 8000 kg d−1, and total riverine discharge varied from 10 to 50 m3 s−1 (Fig. 8). When comparing trends in loading during the different seasons, TP, TN, and discharge were different from each other. There was no temporal trend in TP load and discharge during the spring (Figs. S1.1 and S1.3, Table 1), but there was a declining trend in TN load during spring (Fig. S1.2, Table 1). In summer, TP load and discharge had an increasing trend, while TN load did not show any trend (Figs. S1.4–S1.6, Table 1). Discharge increased in fall (Figs. S1.7–S1.9, Table 1), while in winter, TN and TP load decreased over time (Figs. S1.10–S1.12, Table 1). Examining all three input rates on an annual scale, discharge increased, TP load did not change, and TN load decreased (Figs. S1.13–S1.15, Table 1). For TN, point source load and nonpoint source load both declined statistically significantly. For TP, however, point source load declined but nonpoint source load increased. The Mann-Kendall tests indicate that discharge (τ = 0.165, p = 0.002) was increasing, TN (τ = −0.104, p = 0.062) loads were decreasing, and TP (τ = −0.015, p = 0.504) loads did not change (Fig. 8). Nutrient loading and discharge variables (discharge, TN, TP) were highly correlated (rP = 0.87–0.99; Fig. S6.1) and HVD was moderately correlated with discharge, TN, and TP (rP = 0.46–0.58; Fig. S6.1).

Fig. 8
figure 8

Annual average phosphorus loads (a), nitrogen loads (b), and riverine discharge (c) to the Patuxent River from 1985 to 2021. BFL, below fall line; BFL ps and sho, below fall line point source (ps) + shoreline erosion (sho); and RIM, river input monitoring, those discharges measured directly at the fall line USGS gauging station. The combined dataset is the sum of BFL and RIM

Full size image

Benthic Invertebrate Biomass

We evaluated the impact of HVD on the biomass of benthic invertebrates, hypothesizing that higher HVD would correlate with reduced benthic biomass due to oxygen stress on immobile organisms. At station 71, which is situated in the hypoxic zone of the lower estuary between those two stations with the longest hypoxia duration (LE1.1 and LE1.2; Figs. 1 and 5), we observed a negative relationship between summer benthic biomass and HVD (r2 = 0.34, p < 0.01; Fig. 9). However, no significant relationship was found between spring biomass and HVD. In performing linear regressions between HVD and summer biomass for specific benthic groups, bivalve biomass (r2 = 0.24, p < 0.01) and worm biomass (r2 = 0.14, p = 0.03) exhibited negative relationships with HVD (with p-values below 0.05), suggesting an exponential decrease in biomass with increasing HVD (Fig. 9; Table S2). Conversely, small crustaceans were not substantially related to HVD (r2 = 0.10, p = 0.25, Fig. 9c). Additional linear regressions between discharge and benthic biomass yielded no consistent trends (data not shown).

Fig. 9
figure 9

Relationship between estuary-wide HVD and natural logarithm of total benthic biomass (a), bivalve biomass (b), small crustacean biomass (c), and worm biomass (d) at station 71 in the lower part of the estuary. Each panel illustrates the HVD-biomass relationship, presenting summer data as color-filled circles and spring data as hollow blue circles. Benthic biomass values are represented in natural logarithm scales of grams ash-free dry weight (g DW) on the y-axes. Fitted lines depict the relationships, and corresponding equations, r2, p-values, and sample sizes (n) are provided for each relationship

Full size image

Discussion

Our findings reveal a surprising trend of increasing HVD over the 1985–2021 period, despite a decline in nitrogen loads from the watershed to the estuary. The fact that HVD increased while nitrogen loads decreased (Table 1) contradicts the expectation that reducing nutrient load would lead to spatially constrained or shorter-duration hypoxic volumes. Our results are consistent with the results of Testa et al. (2008), who found higher HVD per unit TN load in the Patuxent over the 1985–2003 period. Variability in hypoxia primarily stems from biological factors, such as water column primary production and respiration (Li et al. 2016), as well as physical factors including river flow and nutrient input (Murphy et al. 2011; Lee et al. 2013; Testa et al. 2022), and both can be modulated by climate change (Najjar et al. 2009; Diaz et al. 2019; Testa et al. 2022). Our analysis of long-term oxygen depletion patterns in the Patuxent River estuary suggests a substantial role of physical factors (river flow, temperature) in overcoming potential improvements associated with nitrogen load reductions. These complex hypoxia dynamics have also been reported for many large estuarine hypoxic zones, including the mainstem of Chesapeake Bay (Lake and Brush 2015; Ni et al. 2020; Frankel et al. 2022), the Baltic Sea (Meier et al. 2018), Long Island Sound (Whitney and Vlahos 2021), and the northern Gulf of Mexico (Laurent et al. 2018). Several prior investigations have suggested that high levels of hypoxia can persist despite nutrient load reductions because hypoxia stimulates nutrient recycling (Testa and Kemp 2012; Savchuk 2018) that can support sustained phytoplankton growth and thus keep a system in a eutrophic regime. We cannot conclude that the Patuxent is in a eutrophic steady state that persists despite nutrient load reductions, in part because the sustained or increasing hypoxia we estimated could simply be a result of long-term changes in external forces (warming, freshwater input) that now control interannual variability.

One possible explanation for increased HVD over time is a long-term increase in freshwater input, which would increase stratification strength to limit bottom water oxygen replenishment in the lower estuary. Nitrogen loads are typically correlated with riverine discharge, which can increase hypoxia (Hagy et al. 2004; Li et al. 2016). In the Gulf of Mexico, models suggest that increasing freshwater input would increase nutrient inputs leading to elevated algal growth and more severe periods of hypoxia (Lehrter et al. 2017). Increased seasonal and annual freshwater flow has been shown to contribute to enhanced stratification within the Patuxent River (Breitburg et al. 2003a). Our trend analysis revealed both long-term increases in freshwater input and fall stratification which could have driven the long-term increase in HVD (Fig. 3). In addition to enhancing stratification, elevated freshwater fluxes, especially during large storms (e.g., 2011, 2018), can deliver large organic matter pulses, which, if retained during winter and subsequently respired in the spring, could initiate hypoxia earlier (Taft et al. 1980; Testa and Kemp 2014). Summer and fall storms may also generate significant wind mixing, episodically disrupting stratification (Breitburg et al. 2003a; Miller et al. 2006a, b) and counteracting the enhanced stratification effect of these flow events. Expectations for hypoxia responses to large storms are therefore complicated by uncertainty in the magnitude of change in these two opposing forces.

Regional climate cycles can also influence hypoxia dynamics through their effects on wind patterns, temperature, and precipitation, as well as the resulting impacts on vertical mixing, phytoplankton productivity, and nutrient loading. For example, past research in the region has associated variability in hypoxic volume in the mainstem Chesapeake Bay with the North Atlantic Oscillation (NAO) through changes in wind direction (e.g., Scully 2010), but large-scale temperature and precipitation patterns in the region within recent decades have not been clearly linked to the El Niño-Southern Oscillation (ENSO) or the NAO (Miller et al. 2006a, b). However, Schulte et al. (2016) did find that ENSO and the Pacific Decadal Oscillation (PDO) were contributing factors to anomalous precipitation variability in the mid-2000s for the Susquehanna River basin, which has some potential impact on the lower Patuxent estuary. In a comparative analysis of the Adriatic and Chesapeake Bay, Brush et al. (2020) reported a strong relationship between chl-a and the NAO, with lower NAO values associated with higher chl-a in the Potomac estuary. Thus, although there is likely to be some degree of influence of these cycles on hypoxia in the Patuxent, trying to discern the magnitude of their effect is outside the scope of our paper, and there is no clear trend in these cycles since 1985.

The observed increases in temperature likely exacerbated hypoxia by reducing the solubility of oxygen in the water column (Najjar et al. 2009) and promoting oxygen-consuming respiratory processes (Testa et al. 2022). For example, the increase of AOU during the fall revealed that there was an elevated biological or physical (e.g., stratification) driver of oxygen declines in the fall across all stations. The elevated AOU could be related to oxygen consumption being enhanced by rising water temperature or by reduced mixing associated with higher summer river flows. Other estuaries worldwide, such as the Baltic Sea (Conley et al. 2009; Carstensen et al. 2014), Danish Straits (Conley et al. 2007), Chesapeake Bay (Ni et al. 2020), the northern Gulf of Mexico (Laurent et al. 2018), and St. Lawrence estuary (Gilbert et al. 2005), have experienced an increase in hypoxia and hypoxic volume, in part due to elevated temperatures that intensify respiration rates. Given that elevated respiration rates lead to increased nutrient regeneration, warming-enhanced respiration could create a feedback loop where nutrient remineralization leads to enhanced phytoplankton productivity, which in turn increases the extent or duration of hypoxia (Najjar et al. 2009). Such feedbacks have been proposed as a result of hypoxia-induced amplification of nitrogen and phosphorus recycling (Conley et al. 2009; Testa and Kemp 2012); thus, the combined effects of warming and elevated hypoxia could increase the likelihood of these feedbacks initiating.

Our estimates of hypoxia onset day suggest that hypoxia is internally generated in the Patuxent River estuary and is thus responsive to local factors. Hypoxia was initiated in the middle of the Patuxent estuary, around station LE1.1, and subsequently occurred at adjacent stations in the lower (LE1.2 and LE1.3) and middle (RET1.1) regions of the estuary. This spatial pattern contrasts with the mainstem of Chesapeake Bay, where hypoxia initiates in northern regions and migrates south as the summer progresses (Testa and Kemp 2014). This region also has the longest duration of hypoxia (Figs. 4 and 5), likely owing to relatively deep waters (13 m) that are strongly stratified during the spring and summer and limit the amount of oxygen that can replenish bottom waters (Boicourt 1992). Pearson correlations between onset day and winter- and spring-averaged bottom water chl-a, stratification strength, and river flow did not explain substantial amounts of variability in onset day across the estuary, which contrasts with studies in the mainstem Bay that found both winter-spring river flow and bottom chl-a strongly co-vary with onset day (Testa and Kemp 2014). We did find a correlation between stratification strength and onset day specifically at LE1.1, suggesting a physical control on hypoxia onset in this region and consistent with the correlation we found between river flow and hypoxic volume days (Fig. S6.1). Future studies could consider the role of other physical factors on hypoxia onset day, such as wind speed and direction (Scully 2010; Wang et al. 2015).

Analysis of estimated water column oxygen depletion rates reveals spatial variability within the estuary related to physical setting. Higher rates of depletion were found in the middle region of the estuary (RET1.1, LE1.1–LE1.4) relative to the tidal fresh region, and these middle stations are those that become hypoxic annually. We explored regressions between depletion rate and winter and spring-averaged bottom water chl-a, stratification strength, and river flow, but these did not explain substantial amounts of variability in depletion rates. This contrasts with studies in the mainstem Bay that found that winter-spring bottom chl-a strongly covaried with depletion rates (Testa and Kemp 2014). When comparing the rate of oxygen depletion in the Patuxent River estuary (0.037–0.083 mg L−1 d−1 or 1.18–2.65 μM d−1) to the mainstem Bay, it falls within the expected range for the season based on previous estimates (Taft et al. 1980; Kemp et al. 1992; Boynton and Kemp 2000; Hagy et al. 2004; Testa and Kemp 2014). The lack of a clear statistical explanation for the Patuxent River estuary’s oxygen depletion rates may be due to conflicting patterns in key controls, most notably an increase in stratification that would enhance depletion rates (Kemp et al. 1992) vs. a decline in TN load that would depress these rates (Fig. 8). Spring and summer chl-a in the middle estuary were stable over time, but highly variable, supporting the notion that the amount of organic fuel for hypoxia did not change.

Finally, the relationships between HVD and benthic biomass in the Patuxent River estuary highlight a crucial ecological connection between oxygen availability and benthic communities. During the summer, hypoxia occurs in the mid and lower estuary, leading to insufficient oxygen levels for benthic organisms to thrive (Breitburg et al. 2003a). In contrast, tidal mixing and low stratification in the tidal fresh estuary maintain adequate oxygen levels (Fig. 2) despite high concentrations of water column nutrients and chl-a (Nezlin et al. 2009, Figs. S5.1–S5.4), conditions that support benthic secondary production throughout the spring and summer. The observed negative relationship between HVD and benthic biomass at station 71 is consistent with prior studies showing reduced macrobenthic production and/or abundance and community structure under hypoxic conditions (e.g., Long and Seitz 2009; Sturdivant et al. 2014). Llansó (1992) also highlighted species-specific tolerance variations and their implications for benthic organism response and recovery from hypoxia events. While our study is currently centered on the broader analysis of three taxonomic groups (bivalves, small crustaceans, and worms), it recognizes the potential for further refinement through a species-specific or functional group approach. The observed influence of HVD on the taxonomic groups highlights distinctions in sensitivity, with bivalves exhibiting the highest response to HVD, followed by worms and then small crustaceans. By delving into specific taxa that dominate biomass, particularly those demonstrating higher sensitivity, future research could uncover nuanced ecological implications. The significance of taxa-specific responses becomes apparent when considering trophic dynamics. Ihde et al. (2015) forage report highlighted benthic/epibenthic invertebrates as comprising six of the ten most important forage taxa by biomass in their study of the mainstem, underscoring the potential cascading effects of benthic biomass reductions on higher trophic levels. Notably, reductions in benthic biomass and abundance expand spatially and temporally with the severity and duration of hypoxic events (Llansó 1992). Our study’s examination of hypoxia onset and breakup adds a temporal dimension to the observed pattern and suggests that the timing of recovery by the benthos following the breakup of hypoxia can vary widely, even in a relatively small system such as the Patuxent River.

Conclusions

Our analysis of long-term changes in metrics of hypoxia and several controlling variables suggests a long-term increase in Patuxent River hypoxia over the 1985–2021 period, which is driven by increasing temperature and freshwater inputs. The fact that HVD increased despite reductions in nitrogen loading and no statistically significant changes in phosphorus loading suggests that mitigating the effects of climate change will require more aggressive nutrient reduction efforts, especially for nonpoint sources (Zhang et al. 2023), to combat increased flow and decreased oxygen concentration associated with climate change. Thus, the concurrent warming and wetting of the Patuxent River estuarine complex present a challenge for management actions targeting eutrophication reduction, and future modeling analysis will help refine our understanding of the balance between eutrophication and climate change in altering the timing, duration, and extend of hypoxia in this estuary. Incorporating climate change into management actions is also critical because of the severe impact of increasing HVD on benthic organisms documented in this study. This supports actions by the Chesapeake Bay nutrient management strategy to seek more aggressive nutrient load reductions in the face of climate change (Cerco 2022). As point source reductions have been documented, nonpoint source reductions especially from the agricultural areas are critical to the restoration of this estuary (Zhang et al. 2023).