
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38266048
202307065
10.1073/pnas.2307065121
early-career-researcherEarly-Career Researcherresearch-articleResearch ArticleecoEcology414
Early-Career-Researcher
Biological Sciences
Ecology
Macroscale controls determine the recovery of river ecosystem productivity following flood disturbances
Lowman Heili E. hlowman@unr.edu
a 1 https://orcid.org/0000-0002-2939-9225

Shriver Robert K. a https://orcid.org/0000-0002-4590-4834

Hall Robert O. Jr. b https://orcid.org/0000-0002-0763-5346

Harvey Judson W. c
Savoy Philip c
Yackulic Charles B. d https://orcid.org/0000-0001-9661-0724

Blaszczak Joanna R. a https://orcid.org/0000-0001-5122-0829

aDepartment of Natural Resources and Environmental Science, University of Nevada Reno, Reno, NV 89557
bDivision of Biological Sciences, Flathead Lake Biological Station, University of Montana, Polson, MT 59860
cU.S. Geological Survey, Earth System Processes Division, Reston, VA 20192
dU.S. Geological Survey, Southwest Biological Science Center, Flagstaff, AZ 86001
1To whom correspondence may be addressed. Email: hlowman@unr.edu.
Edited by Andrea Rinaldo, Ecole Polytechnique Federale de Lausanne, Lausanne, Switzerland; received April 28, 2023; accepted December 4, 2023

24 1 2024
30 1 2024
24 7 2024
121 5 e230706512128 4 2023
04 12 2023
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

River ecosystems rely on varied flows, including regular floods, to provide food and habitat for aquatic organisms. However, flows of freshwater are becoming increasingly managed for irrigation, industry, and other human activities, and the frequency of floods is changing. Our study used time-series data of photosynthesis from 143 rivers across the United States and developed a modeling framework to examine how algae, the base of most riverine food webs, recovered following scouring disturbance during floods. We found that algae in wider rivers recovered more quickly following disturbance (e.g., removal from the bottom) but that the flow thresholds at which algae is disturbed are likely more strongly influenced by site-specific characteristics.

River ecosystem function depends on flow regimes that are increasingly modified by changes in climate, land use, water extraction, and flow regulation. Given the wide range of variation in flow regime modifications and autotrophic communities in rivers, it has been challenging to predict which rivers will be more resilient to flow disturbances. To better understand how river productivity is disturbed by and recovers from high-flow disturbance events, we used a continental-scale dataset of daily gross primary production time series from 143 rivers to estimate growth of autotrophic biomass and ecologically relevant flow disturbance thresholds using a modified population model. We compared biomass recovery rates across hydroclimatic gradients and catchment characteristics to evaluate macroscale controls on ecosystem recovery. Estimated biomass accrual (i.e., recovery) was fastest in wider rivers with less regulated flow regimes and more frequent instances of biomass removal during high flows. Although disturbance flow thresholds routinely fell below the estimated bankfull flood (i.e., the 2-y flood), a direct comparison of disturbance flows estimated by our biomass model and a geomorphic model revealed that biomass disturbance thresholds were usually greater than bed disturbance thresholds. We suggest that primary producers in rivers vary widely in their capacity to recover following flow disturbances, and multiple, interacting macroscale factors control productivity recovery rates, although river width had the strongest overall effect. Biomass disturbance flow thresholds varied as a function of geomorphology, highlighting the need for data such as bed slope and grain size to predict how river ecosystems will respond to changing flow regimes.

ecosystem recovery
disturbance
river
algae
ecosystem metabolism
National Science Foundation (NSF) 100000001 OIA-2019528 Heili E LowmanRobert O. HallJoanna Roberta Blaszczak
==== Body
pmcAnthropogenic change is modifying the frequency and intensity of disturbances in streams and rivers (1). Climate and land use change as well as water extraction and river regulation modify the natural flow regimes of rivers (2) with cascading effects on ecosystem processes and ecological communities (3, 4). As global temperatures rise, extreme climatic events could increase flow variability, including both more regional droughts and more frequent, intense precipitation events (i.e., weather whiplash) (5, 6). Ongoing changes to land use, including urbanization, affect the speed at which runoff is transported into rivers, potentially augmenting or mitigating flood potential (7, 8) and suppressing river ecosystem recovery, specifically autotrophic productivity, following storms (9–11). Water extraction and river regulation (e.g., dams) can further remove or diminish floods, thus homogenizing flows (1, 8) which can alter responses to disturbance and increase autotrophic productivity downstream (12). Few studies have examined how autotrophic productivity in rivers is affected by disturbance at large scales, primarily because these data remain difficult to collect in individual rivers. Therefore, it remains uncertain how variable river recovery may be following high-flow disturbance events. By leveraging increasingly widespread high-frequency, long-term, and large-scale observational datasets (13) and applying process-based models, we can gain insight into controls on river ecosystem recovery at both continental and decadal scales.

In studies of ecosystem disturbance, rates of ecosystem recovery can be quantified from the change in estimates of primary productivity in pre- and postdisturbance states (14, 15). Recovery following disturbance in terrestrial ecosystems is typically measured on a yearly to decadal scale, where disturbance may take the form of drought (16), wildfires (17), warming temperatures (18), and rising seawater levels (19). Due to physiological and metabolic factors as well as more frequent disturbances, aquatic primary producers typically recover faster following disturbance than terrestrial plants, on daily to seasonal timescales [e.g., algae in streams following flash floods (14, 20, 21), giant kelp facing wave disturbance (22)]. To examine changes through time or in the context of a given disturbance, biomass is typically the measured metric of primary production. However, measuring biomass in river ecosystems at high spatial and temporal resolutions is particularly difficult due to its heterogeneous and ephemeral nature, and observations cannot match the spatial extent and frequency for which we have ecosystem metabolism estimates. Instead, estimates of gross primary production (GPP) in rivers (23) can be modeled as a function of underlying autotrophic biomass dynamics (24), enabling the prediction of recovery trajectories in rivers following hydrologic disturbance. Here, we model the biomass of primary producers in rivers, including all photoautotrophs (i.e., periphyton, phytoplankton, and macrophytes), encountering flow disturbance. We use high-frequency estimates of ecosystem primary production and build on prior work illustrating how flow variability relates to riverine GPP (25).

Predictions of postdisturbance recovery of river ecosystem productivity change with the scale of observation and frequency of disturbance. Observations within a river reach suggest increased disturbance frequency may increase recovery rates (20) while whole reach observations show increased disturbance frequency may decrease recovery rates (9). In rivers experiencing relatively high disturbance frequency, predominantly urban or agricultural watersheds, GPP can have low resistance to high-flow events (9, 10). In dammed or otherwise flow-regulated rivers experiencing lower disturbance frequency, management of flow can affect productivity, increasing GPP in steady-flow scenarios and decreasing GPP during periods when flows are altered to meet power demands (12). In addition to factors affecting flow regimes, the relative importance of conditions such as light availability (26) and nutrients (27) may explain variation in recovery rates because these conditions facilitate high rates of algal productivity when flows are steady. River width tends to increase with the size of the drainage area (28). As a result, less of the surface area tends to be shaded by tree canopy, increasing light availability (29). However, an increase in light availability may not translate to greater light penetration through the water column because many wider rivers are deep and often have high turbidity (30). River size may also modulate the effects of land use; in particular, wide rivers dampen the effects of watershed development on productivity dynamics (31). Nutrient availability in rivers is not thought to be a primary driver of variation in GPP (26), but river productivity may increase in human-altered watersheds due to alleviation of nutrient limitation (31). Each of these interactions makes it difficult to predict how factors such as river size and nutrient availability may determine the rate of recovery of autotrophic productivity.

In this study, we employ a macrosystems approach (32) to examine the relative influence of cross-scale controls (i.e., from the river reach to landscape scale), in addition to light and flow, that mediate the relationship between disturbance frequency and recovery rates of autotrophic productivity following flood disturbances in rivers across the continental United States. We modify a state-space time series modeling approach developed by Blaszczak et al. (24) that uses daily GPP estimates (23) and a density-dependent population model (33) to predict daily autotrophic biomass and estimate recovery following high-flow disturbance events. We then investigate macroscale controls of ecosystem recovery, including land cover, river size, temperature, and nutrients, following flood disturbance across varied hydrologic regimes. In addition, we estimate the magnitude of flow which causes reductions in productivity through time (i.e., the disturbance flow threshold) and compare river-specific thresholds to bed disturbance thresholds based on the geomorphic characteristics of rivers. Ultimately, we quantify the recovery of river productivity following flood disturbances to better inform predictions of how river ecosystems spanning a wide range of hydrologic and geographic conditions will respond to changing disturbance regimes.

Results

Estimates of Maximum Growth Rates and Biomass Accrual.

Model-estimated maximum growth rates ( rmax ) of autotrophic biomass were relatively fast but varied widely across all rivers (n = 143). Most rivers had a doubling time averaging 1 wk (7.3 d doubling time; median rmax across all sites = 0.095 d−1), with few rivers having exceptionally fast growth rates resulting in a right-skewed rmax estimate distribution. The eight sites that fell at or above the 95th percentile rmax value (0.35 d−1 or 2.0 d doubling time) displayed lower coefficients of variation in discharge (mean of eight sites = 0.86), wider river size (mean width across eight sites = 64 m), and higher daily GPP (mean of eight sites = 9.4 g O2 m−2 d−1). Similarly, the eight sites whose estimated rmax values fell at or below the 5th percentile (0.028 d−1 or 24 d doubling time) instead displayed higher coefficients of variation in discharge (mean of eight sites = 1.74), narrower river size (mean width across eight sites = 10 m), and lower daily GPP (mean of eight sites = 0.5 g O2 m−2 d−1). Overall, uncertainty intervals of rmax estimates were larger at higher values (maximum rmax estimate across all sites = 0.66 d−1 or 1 d doubling time, CI95 = [0.54, 0.79]) and smaller at lower values (minimum rmax estimate across all sites = 0.006 d−1 or 124 d doubling time, CI95 = [−0.044, 0.052]).

The recovery rates of river biomass following disturbance, as estimated by maximum biomass accrual rates ( amax   ), covaried more with in-stream conditions and less with flow modifications. Similar to rmax   estimates, few rivers had very fast recovery rates following flood disturbance, resulting in a right-skewed amax   distribution. Higher values (i.e., above the 95th percentile) of amax   also corresponded to lower coefficients of variation in discharge, wider river size, and higher daily GPP while lower amax   values (i.e., below the 5th percentile) corresponded to higher coefficients of variation in discharge, narrower river size, and lower daily GPP (SI Appendix, Fig. S6). Multilevel models constructed for amax   values showed that river width, and by extension watershed size, had the strongest effect on biomass accrual (Fig. 1 and SI Appendix, Fig. S5, n = 137). amax increased, but in an attenuating way with river width ( βlog10(Width) = 0.47, CI95 = [0.27, 0.67]). A site with more annual exceedances of the estimated critical discharge disturbance threshold ( Qc ) also predicted higher amax values ( βExceedances = 0.18, CI95 = [0.03, 0.33]). The certain presence of an upstream dam predicted lower amax values ( βDam = −0.40, CI95 = [−0.71, −0.09]). Neither increasing water temperature nor road density in the watershed strongly affected variation in amax values, nor did including dissolved nutrient concentrations in the second model structure ( βTemperature = 0.00, CI95 = [−0.17, 0.19]; βRoad = 0.01, CI95 = [−0.19, 0.22]; βlog10(NO3) = 0.16, CI95 = [−0.11, 0.43]; βlog10(P) = 0.12, CI95 = [−0.16, 0.40]).

Fig. 1. (A) Multilevel regression model estimates indicate that three covariates—dams, disturbance threshold exceedances, and river width—strongly covaried with maximum river biomass accrual ( amax   ). Covariates are colored according to the information they provide (in-stream conditions in blue and flow modifications or regimes in orange). Median parameter estimates are denoted by points, and 95% credible intervals are denoted by shaded lines. These results suggest that accrual (B) decreased with increasing interference by dams upstream of a site, had no clear relationship with (C) mean water temperature or (D) road density in a watershed, (E) increased exponentially with increasing mean annual disturbance threshold ( Qc   ) exceedances, and (F) increased in an attenuated manner with increasing river width. In panels B–F, shaded points show individual site data (n = 137), and if a covariate had a strong effect on accrual (i.e., 95% CIs did not cross 0), panels also present 100 draws of the joint distribution of the linear model.

Estimates of Flow Disturbance Thresholds.

Flows that scoured biomass were widely variable across all rivers (n = 130). The distribution of critical disturbance flow thresholds ( c ) was uniformly distributed, with a median value of 0.55 which indicates that the critical disturbance flow threshold of biomass is 55% of the magnitude of maximum discharge measured during a given time series. Values of c also varied 100-fold across the 130 sites. The seven sites that fell at or above the 95th percentile c value (1.08 or 108% of the magnitude of maximum discharge) displayed much lower coefficients of variation in discharge (mean of seven sites = 1.65), whereas sites whose estimated c value fell at or below the 5th percentile (0.043 or 4.3% of the magnitude of maximum discharge) instead had higher coefficients of variation in discharge (mean of seven sites = 3.87). Neither river width (95th percentile mean width = 28 m; 5th percentile mean width = 11 m) nor daily GPP (95th percentile mean = 2.34 g O2 m−2 d−1; 5th percentile mean = 1.87 g O2 m−2 d−1) covaried meaningfully with c values as rmax or amax values did. Unlike rmax estimates, uncertainty intervals of c parameter estimates did not have as clear of a relationship with the magnitude of the estimate.

Biomass disturbance flows ( Qc   ) were mostly lower than the 2-y flood ( Q2yr   ). Seventy-three percent of sites (n = 95) had Qc:Q2yr   < 1, indicating that their biomass disturbance threshold fell below the 2-y bankfull flood, although several rivers had Qc:Q2yr   ≫ 1 (maximum estimated Qc:Q2yr   = 6) creating a right-skewed Qc:Q2yr   distribution. Unlike c   values, Qc:Q2yr   values did not covary with the coefficient of variation in discharge. While higher Qc:Q2yr   values were observed in rivers that were wider and less productive, the differences between high and low Qc:Q2yr   value sites were less exaggerated than the trends displayed with c   values (SI Appendix, Fig. S7). None of the three predictors investigated in the multilevel model—river width ( βlog10(Width)   = 0.01, CI95 = [−0.22, 0.25]), road density ( βRoad   = −0.07, CI95 = [−0.30, 0.17]), or interference by dams ( βDam   = −0.19, CI95 = [−0.57, 0.20])—had a meaningful covariance with Qc:Q2yr   values, indicating that river width, land use, and watershed characteristics did not strongly affect biomass disturbance thresholds (Fig. 2, n = 124).

Fig. 2. (A) Multilevel regression model estimates indicate that none of the covariates investigated have a strong effect on biomass disturbance thresholds scaled to a site’s 2-y flood ( Qc:Q2yr ). Covariates are colored according to the information they provide (in-stream conditions in blue and flow modifications or regimes in orange). Median parameter estimates are denoted by points, and 95% credible intervals are denoted by shaded lines. Values of Qc:Q2yr show little change with differences in (B) likelihood of interference by dams, (C) road density in a watershed, or (D) river width. Panels B–D present shaded points showing individual site data (n = 124).

Comparison of Model-Estimated and Geomorphically Estimated Flow Disturbance Thresholds.

Our model-estimated Qc:Q2yr values (i.e., normalized biomass disturbance thresholds) were greater than all of the normalized bed disturbance thresholds ( Qg:Q2yr ) estimated using the geomorphic model, except at one site, Accotink Creek, which had the steepest slope and the largest median grain size (Table 1). The biomass disturbance thresholds at two sites, Kankakee River and Difficult Run, were 7 and 34 times larger than the bed disturbance thresholds at those sites, respectively. In particular, the biomass disturbance threshold at Difficult Run was the greatest, but the bed disturbance threshold at this site was the least out of all five sites. Of the three wider rivers examined (i.e., South Fork Iowa River, Iroquois River, and Kankakee River), all three displayed greater biomass than bed disturbance thresholds, and the relative difference between these values increased with decreasing median grain size and bed slope.

Table 1. Model estimated autotrophic biomass disturbance thresholds (Qc ) were typically greater than bed disturbance thresholds (Qg) once normalized to a site’s 2-y flood (Q2yr)

Site	Median grain size (mm)	Bed slope (m m−1)

	Site description	Biomass model	Geomorphic model	
Biomass disturbance threshold, Qc (m3s−1)	Biomass disturbance threshold, Qc/Q2yr	Bed sediment movement threshold, Qg (m3s−1)	Bed disturbance threshold Qg/Q2yr	
Accotink Creek, VA* (GPPdaily = 0.65 g O2 m−2 d−1)	16	0.005	4-m wide, pool and riffle, incised with restoration cross-vanes, partly shaded	1.6 [1.08, 5.45]	0.44	1.9	0.51	
S. Fork Iowa River, IA (GPPdaily = 2.66 g O2 m−2 d−1)	11	0.001	20-m wide riffle and run, gravel-bed, mostly unshaded	41.6 [39.4, 50.3]	0.84	28.8	0.58	
Difficult Run, VA (GPPdaily = 0.45 g O2 m−2 d−1)	6	0.005	3-m wide, pool and riffle, gravel-bed, shaded	14.5 [2.1, 18.2]	2.37	0.4	0.07	
Iroquois River, IN (GPPdaily = 0.68 g O2 m−2 d−1)	2	0.00018	20-m wide riffle and run, grave/sand-bed, partly shaded	69.7 [60.7, 72.7]	0.77	13.5	0.15	
Kankakee River, IN (GPPdaily = 0.50 g O2 m−2 d−1)	0.4	0.00037	20-m wide dune and ripple, sand-bed, partly shaded	50.8 [40.9, 64.0]	1.23	7.5	0.18	
Beneath each site name is its mean daily gross primary production (GPP), followed in the next columns by median grain size, bed slope, and a general site description.Qcvalues are reported with 95% CIs in brackets underneath. Note that sites denoted with an asterisk passed our model convergence diagnostics ( R∧< 1.05) for maximum growth rate estimates ( rmax ) but not for estimates of critical discharge thresholds at which autotrophic biomass is disturbed ( c).

Discussion

Several macroscale controls were strong predictors of river productivity recovery as represented by biomass accrual rates ( amax ) in river sites across the United States. Narrower rivers with less frequent disturbances and a higher likelihood of upstream dams had slower recovery of river productivity following flow disturbances. We did not find evidence of similar macroscale controls on biomass disturbance thresholds ( Qc:Q2yr ), but estimated Qc:Q2yr values were higher than geomorphically derived bed disturbance thresholds ( Qg:Q2yr ). This finding suggests that the biomass model is more sensitive to large disturbances of river primary productivity with higher ecosystem disturbance threshold discharge values compared to the geomorphic disturbance threshold discharge values. In other words, the geomorphic threshold for bed disturbance may underestimate the threshold for large disturbances to river primary productivity and may instead estimate smaller, shorter-lived disturbances. Together, these findings indicate that while larger disturbances and longer recoveries of river productivity may be controlled by multiple, interacting macroscale factors, smaller disturbance events may depend more specifically on geomorphology including factors such as bed slope and grain size.

Biomass Accrual and Implications for Recovery Following Disturbance.

Recovery of river productivity, as quantified by maximum biomass accrual or amax   , varied as a function of both river and watershed characteristics. Across 137 river reaches, wider rivers had higher amax   values and recovered more quickly following flow disturbance. Wider, higher order rivers have greater light availability and steadier flows, conditions that maximize annual productivity (26) and, according to our results, recovery of the autotrophic biomass that underlies ecosystem productivity. In addition to in-stream conditions, watershed characteristics and metrics of flow variability covaried with river recovery. amax   increased in sites with lower likelihood of dams upstream. Dams may dampen seasonal variability and homogenize flows (34), and while tailwaters can have high rates of primary production, they can shift the primary location of production from benthic to planktonic, due to productivity in upstream reservoirs (35). Flows beneath dams may vary at fine temporal scales due to closely following daily electrical power demand (i.e., “load-following flows”); in these circumstances, with high daily flow variability, productivity decreases due in part to increases in turbidity (12). At our sites, the presence of upstream dams negatively related with amax   and suppressed recovery of river productivity, suggesting that the flow beneath these dams was perhaps more variable and less conducive to higher productivity conditions than larger dams, such as those on the Klamath and Colorado Rivers (12, 35). amax also increased with the number of annual high-flow events exceeding the biomass disturbance threshold, which suggests that autotrophic populations recovered more rapidly in more frequently disturbed rivers. Although rivers with higher coefficients of variation in discharge displayed slower recovery rates, these locations were also narrower, less productive rivers, and our model suggests that river width had a stronger effect on overall recovery than biomass disturbance threshold exceedances alone. River algal assemblages can depend on flow disturbance severity (36), thus more frequently disturbed rivers may select for rapidly growing, disturbance-adapted algal species. Of all river and watershed predictors with effects on recovery, both river width and dams upstream have demonstrated relationships with flow variability, suggesting that changes in flow regime likely control the rate at which river productivity recovers. However, we suggest that river characteristics, namely river width, most strongly affected the recovery of autotrophic productivity following high-flow disturbance events.

Water temperature, watershed land use, and nutrient concentrations were not related to recovery rates following disturbance events. Few studies have documented a strong link between temperature and riverine productivity (37), and the direct effect of temperature on productivity recovery times may be difficult to estimate due to the effect of temperature on other in-stream processes that may indirectly control primary productivity, such as grazing (38), light, or nutrient availability (39). In particular, grazing by algivores is not accounted for in this modeling framework, although consistent differences among sites would implicitly be included in rmax   and carrying capacity ( K   ) estimates. Scouring flood events may alter consumer (i.e., fish) mediated effects on algal growth (40, 41), and these internal processes may explain some of the variance unexplained by our models; however, these data are not available at a macroscale. Increasing developed land cover (i.e., road density) may lower river metabolism due to less infiltration of precipitation and consequently more frequent, flashier flows that scour streambeds and remove biomass (9, 11). Watershed road density did not covary with amax   , but variability in amax   estimates was high at lower road densities, suggesting that other attributes of watershed development and land use may influence autotrophic population dynamics. Because nutrient data were only available at a third of the sites (n = 50), we are unable to comment further on the lack of relationship between nutrients and amax   except to say that water column nutrient concentrations did not appear to control recovery rates. Other studies have also found that river nutrient concentrations do not predict nutrient limitation (42, 43), which suggests nutrient concentrations alone may not explain variation in GPP. Overall, we propose that in-stream, physical attributes (i.e., river width) that affect flow and bed stability have a stronger influence than abiotic conditions (i.e., temperature, nutrients) on how rapidly rivers recover following high-flow disturbance events.

Aquatic ecosystems have more frequent disturbances (e.g., stormflows) and primary producers (i.e., macrophytes, periphyton, phytoplankton) with much shorter lifespans than in terrestrial ecosystems, reducing the timescales (hours to seasons) at which the ecosystem can recover from disturbances (44, 45). Across multiple studies in streams and rivers, the mean recovery time of GPP was 7 d (10, 11, 46–48), which matches our median estimated biomass doubling time (i.e., rmax ) of 7.3 d and suggests the time periods of recovery estimated by our model structure were biologically reasonable. While aquatic algal maximum growth rates and the corresponding doubling times may be on the order of days to weeks, the doubling time of terrestrial plants such as sedges, herbs, shrubs, and trees instead ranges from weeks to decades (49). The few sites (n = 8) whose estimated recovery times were longest, weeks to months, may be more representative of macrophyte-dominated rivers or rivers in which algae may colonize the surface of slower-growing macrophytes (50). Still, these rivers likely did not approach longer terrestrial recovery timescales because algal cells can grow more quickly due to their lower tissue thickness (49) and general lack of rigid structural compounds, like lignin (51). However, these rivers did display 10-fold longer recovery times as well as high coefficients of variation in discharge, so recovery of river productivity may have been suppressed by high turbidity (29), which our model was unable to account for and for which data is generally more sparse.

Biomass Disturbance Thresholds and Indications of Disturbance Frequency.

None of the covariates included in our models were strongly related to variability in biomass disturbance thresholds ( Qc   ), even when normalized to the 2-y flood ( Q2yr   ) in each river. We chose Q2yr   as the metric of flow comparison between sites for multiple reasons. First, Q2yr   may represent a coarse estimate of the flow magnitude that mobilizes the most sediment over time (52). In addition, Q2yr   typically represents the recurrence interval of a bankfull discharge, or a flow at which the river reaches but does not overflow its banks (53, 54). Previous analyses of biomass disturbance thresholds also found no relationship between Qc:Q2yr and measures of river size (i.e., watershed area, stream order) (24). Furthermore, neither urban development within a watershed (i.e., road density) nor flow regulation (i.e., dams) covaried with biomass disturbance thresholds. However, Qc remained below Q2yr at 73% of sites (n = 95 out of 130), corroborating a previous study in which all Qc estimates fell below Q2yr values (n = 6 total) (24). Together, these findings indicate that floods regularly scour autotrophic biomass in streams, suggesting generally low resistance to high-flow events.

To provide additional context to our estimates of biomass disturbance thresholds ( Qc   ), we compared our model estimates of Qc   to bed disturbance thresholds ( Qg   ) in five rivers and found that higher discharge was required to disturb biomass than to initiate disturbance of the river bed (Table 1). This finding makes sense if the lower geomorphic flow thresholds, which estimate threshold of movement of median grain size, do not greatly disturb autotrophic biomass. Rather, our modeling outcomes suggest that more than the initiation of bed movement is typically needed to actively scour and substantially remove autotrophic biomass to reset GPP. In limited direct comparison, only Accotink Creek had Qg   > Qc   , but it was also the only highly urbanized site, and it had large erosive peakflows, bed incision, and bed coarsening. Accotink Creek has increased light because of riparian tree removal to facilitate construction of restoration cross vanes designed to lower shear stress along the banks (55) with algae attached to large cobbles that are substantially larger than the median grain size (56). Difficult Run had the highest Qc   relative to the 2-y flood, which suggests that the autotrophic assemblage at that site was particularly resistant to high-flow disturbance in part because of substantial benthic periphyton in channel side-storage zones that were protected from highest shear stresses (56). In the wider rivers, including the South Fork of the Iowa River, the Iroquois River, and the Kankakee River, sites with decreasing median grain size and bed slopes had increasingly larger differences between Qc   and Qg   values, which may indicate the limitations of the geomorphic disturbance threshold as a predictor of substantial disturbance of biomass in wider, finer-grained rivers. Our model is unable to distinguish between high turbidity events and scouring flows, both of which may decrease GPP (9, 29, 57), but our comparison of disturbance thresholds at these five sites highlights the need for data on physical river characteristics, namely slope and grain size, because they may determine disturbance thresholds for riverine autotrophic productivity. More widespread and publicly available characterization of river geomorphological characteristics would improve our ability to predict at continental scales how river ecosystems will respond to changing flow regimes.

Conclusions

Estimates of river metabolism and recovery incorporate both internal ecological dynamics and external forcing from the contributing watershed and enable predicting ecosystem response to disturbance (24). Our estimated metrics of biomass flow disturbance thresholds and river recovery confirm that substantial disturbance occurs frequently in rivers and the biotic response to these disturbance events depends on both river characteristics, such as width and grain size, as well as watershed characteristics, such as upstream dams. In addition, our examination of biomass disturbance and recovery in 143 rivers demonstrates that disturbance thresholds and recovery rates were highly variable. Others have advocated for the study of controls on autotrophic biomass and stream metabolism due to potentially nonlinear relationships with abiotic variables (58). We suggest that additional research should clarify the effects of controls such as land use (e.g., ref. 59), nutrient addition (e.g., ref. 60), and temperature fluctuation (e.g., ref. 61) on river GPP dynamics, the results of which could be used to inform macroscale studies such as ours. Future work focused on modeling river recovery could assess additional controls on river productivity, including turbidity and grazing, and identify site types at which disturbance signals may be challenging to detect, such as spring-fed streams. Furthermore, this study highlights how the effects of changing flow regimes on river ecosystem recovery may depend on the geomorphic context in which disturbances occur, since in-stream conditions, namely river width, were the strongest predictors of recovery rates. As human activity and development increasingly alters river flow and disturbance regimes, it will be necessary to develop approaches that quantify recovery of river productivity integrating both natural and human influences so that we may better predict macroscale patterns and site-specific changes in river ecosystem function.

Materials and Methods

Data Assembly.

We sourced daily estimates of river gross primary productivity (GPP; g O2 m−2 d−1) and associated daily covariates (i.e., discharge) from a published dataset of ecosystem metabolism in 356 rivers between 2007 and 2017 (23). Prior to fitting our state-space time series models, we filtered the Appling et al. (23) GPP dataset to ensure that our inferences were based on high-quality data. We established filters based on model convergence, biologically realistic GPP estimates, site type (e.g., stream, canal), and data continuity; additional details regarding filtering thresholds may be found in SI Appendix.

On days for which we had GPP estimates, we relativized all mean daily discharge values at a site to the maximum observed mean daily discharge at that site, so that all discharge values fall between 0 and 1 which improves model convergence. We chose to normalize discharge to the maximum observed discharge rather than the 10-y flood because model estimates were the same and displayed less variation (SI Appendix, Fig. S8). We paired these discharge time series with an existing dataset of cumulative daily light availability at the stream surface (mol m−2 d−1; 62) which takes into account cloud and canopy cover, and we relativized all light values to the maximum light measurement at a river, so all values fall between 0 and 1. Light estimates were unable to be calculated at several sites (n = 17) due to lack of available azimuth or leaf area index data, and therefore these sites were removed from the dataset. Once the available discharge and light datasets were combined with the filtered daily GPP estimates, the dataset consisted of 181 sites or 684 site-years of data, with daily time series ranging from 3 mo to 9 y in length (median: 870 d), that were fit using the model structure described below (SI Appendix, Figs. S1 and S2 and Table S3). We felt confident using a range of time series lengths because midsize flood events, rather than more rare, extreme flood events that may require longer time series to capture, are most likely to have the greatest cumulative effect on ecological function in rivers (3).

To compare our model estimated parameters of biomass growth and disturbance thresholds with site-specific information, we paired each site with a dataset of site characteristics created by Blaszczak et al. (63), largely sourced from the National Hydrography Dataset Plus V2 (64) and the StreamCat database (65). We accessed nutrient concentration data measured by the U.S. Geological Survey (USGS; 66) using the dataRetrieval R package (67) and filtered between 2007 and 2017, to match our GPP data time frame, for records of dissolved nitrate (USGS parameter codes 00618 and 71851) and phosphorus (USGS parameter codes 00666 and 01072).

To further evaluate model estimates of disturbance thresholds, we compared our estimated disturbance flow thresholds with the 2-y flood ( Qc:Q2yr ) to examine how sensitive autotrophic biomass was to disturbance relative to a rough estimate of mobilization of the entire bed ( Q2yr ). Since the discharge data published by Appling et al. (23) spanned 2007 to 2017, we compiled mean daily discharge data (USGS parameter code 00060) between 1970 and 2020 for each site from data measured by the U.S. Geological Survey (66) using the dataRetrieval R package (67). After removing provisional data, we calculated Q2yr at each site by fitting a linear regression between log maximum annual discharge and the exceedance probability and used the slope of the relationship to calculate the 2-y flood recurrence interval discharge (68). Note that log transformation in this text refers to the natural logarithm ( loge ) unless otherwise stated.

Although detailed grain size data do not exist for all sites included in this study, we could identify five streams and rivers where sufficient grain size data were available. For these five sites, we calculated a threshold for bed disturbance (i.e., the critical Shields stress criterion for initiation of bed movement) to compare how sensitive autotrophic biomass is to disturbance ( Qc   ) relative to bed disturbance determined by geomorphic measures ( Qg   ) (69, 70). While Q2yr is assumed to represent mobilization of the entire bed, our calculation of Qg instead represents the initiation of bed movement. Additional details regarding these calculations are provided in SI Appendix. Model results compared to both thresholds for bed and biomass disturbance in the South Fork Iowa River are displayed in SI Appendix, Fig. S9.

Model Fit.

Following data aggregation, we fit a latent biomass time-series model from Blaszczak et al. (24) to the dataset of daily GPP, discharge, and light. The model structure combines a state-space time series framework with a density-dependent population model, an approach that allows estimating ecologically relevant parameters, such as maximum growth rate, and their uncertainties (71, 72).[1] gt=LteBt+εobs,t

The observation model predicts GPP ( g ) at a daily time step ( t ) as a function of light availability at the stream surface ( L ) and the log daily autotrophic biomass ( B ) while also estimating observation error ( εobs ), which is normal (0,σobs) (Eq. 1). Although aquatic GPP is not strictly linearly related to biomass (73), we chose not to adapt a more complex form so as to apply the model structure to a wide number of rivers and a variety of autotrophic assemblages.[2] Bt=Bt-1+rmax+λeBt-1Pt+εproc,t

The process model predicts the log biomass ( B ) using a density-dependent population model (33). The modified Ricker model structure predicts B as a function of the log of the previous day’s biomass ( Bt-1 ), the maximum growth rate ( rmax ), lambda ( λ=-rmaxK where K is the carrying capacity), and an estimated daily persistence term ( P , Eq. 3) in addition to estimating process error ( εproc ), which is normal (0,σproc) (Eq. 2).[3] Pt=e-e100s(Qt-c)

The second part of the process model predicts the daily persistence term ( P ) as a function of discharge ( Q ), the sensitivity of the persistence transition from presence to removal of biomass ( s ), and the critical disturbance flow threshold at which autotrophic biomass is removed ( c ) (Eq. 3).

To account for larger gaps in time series (i.e., a new sequence of observations beginning in the summer of the following year), we developed a reinitialization process for the process model, the details of which are provided in SI Appendix. For additional information regarding the modeling approach, particularly the development of the process equations and parameter priors as well as original training and testing datasets used, we direct readers to the manuscript and Supporting Information provided in Blaszczak et al. (24).

We fit the model at each individual site (n = 181) using Bayesian inference. Although we investigated fitting a multilevel model to pool information across sites, we ultimately found this approach to be too computationally intensive and the run time prohibitive for enabling a reproducible workflow (approximately 1 to 2 mo to fit the model at all sites on a high-performance computing cluster) and thus fit each site independently. We assigned all parameters slightly informative priors to limit posterior results to realistic values, e.g., positive maximum growth rates (initial value and prior assignments in SI Appendix, Table S1). These values also encouraged model convergence at sites where flow disturbances were few or low in disturbance magnitude, and we based these values on the results of the model development by Blaszczak et al. (24).

We fit latent biomass models with the STAN (74) probabilistic programming language within R (v 4.1.2) (75) using the rstan package (76) in RStudio (v 2021.9.1.372) (77). We ran three chains per site run for 5,000 iterations including 2,500 warm-up or burn-in iterations (i.e., 7,500 total iterations included in parameter summaries presented below) and a maximum tree depth of 12. We fit all models using the Beartooth high-performance computing cluster housed at the University of Wyoming (78). Additional R packages used to extract, visualize, and analyze data include the tidyverse (79), here (80), shinystan (81), ggbreak (82), and patchwork (83) packages. Model results for 16 sites of varying time series length, land cover, light, and flow conditions are displayed in SI Appendix, Figs. S3 and S4.

Model Diagnostics.

Following model fitting, we used model diagnostics and seasonal variation in time series to filter parameter estimates prior to post hoc analyses. First, we removed sites which had potential scale reduction statistic ( R∧ ) values >1.05 for the rmax parameter, negative median rmax estimates, and time series spanning <6 mo (n = 38), resulting in 143 sites remaining that could be included in the regression analysis examining biomass recovery as the dependent variable (SI Appendix, Table S3). We performed an exponential regression of c versus s parameter estimates for each site and found that 18% of sites (n = 32 of 181) had an R2 value >0.4, indicating that these sites had high covariance between c and s parameter estimates (SI Appendix, Fig. S10). To remove sites at which c and s parameters might covary and models that displayed poor convergence, we used the rmax filtered dataset and removed sites at which the c parameter R∧ values were >1.05 (n = 13), resulting in 130 sites with parameter estimates that could be included in the regression analysis examining biomass disturbance thresholds as the dependent variable (SI Appendix, Table S3).

Results Reporting and Site Trends.

Maximum growth rates are reported as both daily estimates (d−1) as well as doubling times (days), which were calculated as log(2)rmax . To better compare results across sites, we used rmax to estimate maximum biomass accrual ( amax ). Maximum accrual is defined as the maximum rate at which a population, in this instance the autotrophic population in the river reach contributing to the dissolved oxygen signal used to estimate ecosystem productivity, may increase without density limitation. Greater accrual estimates indicate greater population growth meaning faster recovery following flow disturbances. We calculated amax as a metric of recovery because it 1) accounts for correlation between model estimates of rmax and λ , 2) is a population-level metric as opposed to describing an individual unit of biomass (i.e., rmax ), and 3) is an integrative metric of postdisturbance recovery. To quantify amax at each site, we modified an equation quantifying maximum sustainable yield for a Ricker model initially developed by Hilborn (84) (Eq. 4).[4] amax=rmax(0.5-0.07rmax)-λ

To compare flows necessary to disturb biomass across sites, we used the model-estimated parameter of the critical disturbance flow threshold ( c ) to estimate the biomass disturbance threshold ( Qc ) relative to the bed mobilization threshold, or the 2-y flood ( Q2yr ). We converted estimated critical disturbance flow thresholds ( c ) to discharge values (m3 s−1) using the maximum discharge ( Qmax ) by which discharge data had previously been normalized. Then, we divided the biomass disturbance threshold ( Qc ) by the 2-y flood ( Q2yr ) at each site, so that values would be comparable across sites. This procedure provides a more meaningful comparison of flow disturbance because magnitude of discharge scales with river size and stream order. Ratio values >1 indicate that the biomass disturbance threshold is greater, and likely occurs less frequently, than the 2-y flood, whereas values <1 indicate that the biomass disturbance threshold is less, and more frequent, than the 2-y flood.

To examine relationships between estimated biomass accrual and disturbance threshold parameters ( amax , Qc:Q2yr ) and site-specific characteristics, we constructed a series of Bayesian generalized multilevel models. A list of all covariates initially considered in all model structures is available in SI Appendix, Table S2 in addition to SI Appendix, Figs. S6 and S7. We initially screened all potential covariates for cross-correlation prior to inclusion in any of the models and log10-transformed amax values, Qc:Q2yr values, and necessary covariates due to left or right-skewness. We also scale-transformed all continuous, numerical covariates prior to their inclusion in the models. For amax values, final model structure included mean water temperature, road density, the likelihood of interference by dams, the mean annual times Qc was exceeded, and river width as completely pooled effects to investigate the effects of climate, surrounding land use, flow regulation, flow regime, and river size, respectively, on biomass accrual (n = 137, SI Appendix, Table S3). The likelihood of interference by dams derived from Appling et al. (23) who evaluated watersheds for metabolism modeling suitability; the 0th, 50th, 80th, and 95th percentiles of the mean 80% oxygen turnover distance were compared to the distance to an upstream dam feature, and the corresponding percentile distance is reported. We combined these into two categories (5 to 50% and 100%), with the 0th percentile indicating the closest distance or a 100% likelihood of interference by a dam. A group-level intercept of the hydrologic unit code denoting a site’s geographic region (HUC2) as defined by the USGS was included to account for the potential effects of a shared river network. An alternative model structure including dissolved nitrate and phosphorus as two additional completely pooled effects was also created to investigate the effect of nutrient availability, but this model was evaluated separately due to limited data availability (n = 50, SI Appendix, Table S3). For Qc:Q2yr values, the model structure included road density, the likelihood of interference by dams, and river width as completely pooled effects to investigate the effects of surrounding land use, flow regulation, and river size on disturbance thresholds (n = 124, SI Appendix, Table S3). A group-level intercept of the HUC2 site region was again included to account for the potential effects of a shared river network.

We fit the model in a Bayesian framework using the brms package (85), with four chains run for 2,000 iterations including 1,000 warm-up iterations (i.e., 4,000 total iterations per model fit) and a maximum tree depth of 12. We assessed model convergence by first visually inspecting for good mixing of chains and lack of divergent transitions as well as ensuring R∧ < 1.05 and the effective sample size >10% for all estimated parameters. We validated all versions of the model structure for best model fit using leave-one-out cross-validation (86) alongside ecological information added by each covariate. Final model results, including median posterior values and 95% credible intervals, were extracted and visualized using the tidybayes package (87) and are displayed in Figs. 1 and 2 and SI Appendix, Fig. S5. We investigated how much the length of the time series included in these analyses influenced our model results and found no relationship between time series length and the residuals of the regression models (SI Appendix, Fig. S11) or parameter estimates (SI Appendix, Fig. S12). We further investigated including only sites with equivalent data density, but found the small number of such sites available (n = 24) provided insufficient power with which to perform a macroscale analysis (SI Appendix, Figs. S13 and S14).

Supplementary Material

Appendix 01 (PDF)

Click here for additional data file.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government. We sincerely thank Ana Miller-ter Kuile, Kiona Ogle, and members of the Modelscape Consortium for their assistance with model development, refinement, and interpretation as well as Alex Buerkle and the entire support staff of the Advanced Research Computing Center's High Performance Computing Cluster at the University of Wyoming, on which all models were run. In addition, we thank Mike Vlah for data support and Leon Katona, Nick Marzolf, Kelly Loria, Jordan Zabrecky, and members of the Blaszczak lab for project narrative refinement. This manuscript was greatly improved by the comments of Erin Hotchkiss, two anonymous reviewers, and the handling editor. Financial support for this project was provided by the Modelscape Consortium, funded by the NSF’s Established Program to Stimulate Competitive Research award #OIA-2019528.

Author contributions

H.E.L. and J.R.B. designed research; H.E.L. performed research; H.E.L., R.K.S., R.O.H., J.W.H., P.S., C.B.Y., and J.R.B. contributed new reagents/analytic tools; H.E.L. analyzed data; and H.E.L., R.K.S., R.O.H., J.W.H., P.S., C.B.Y., and J.R.B. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

Model estimates and all R and STAN scripts, including links to download original datasets, used in the analyses described above are publicly available at https://zenodo.org/doi/10.5281/zenodo.10455822. GPP data were archived by Appling et al. (23). Previously published data were used for this work (23, 29, 62, 63, 66).

Supporting Information

This article is a PNAS Direct Submission.
==== Refs
1 N. L. Poff, B. P. Bledsoe, C. O. Cuhaciyan, Hydrologic variation with land use across the contiguous United States: Geomorphic and ecological consequences for stream ecosystems. Geomorphology 79 , 264–285 (2006).
2 K. Eng, D. M. Wolock, D. M. Carlisle, River flow changes related to land and water management practices across the conterminous United States. Sci. Total Environ. 463–464 , 414–422 (2013).
3 M. W. Doyle, E. H. Stanley, D. L. Strayer, R. B. Jacobson, J. C. Schmidt, Effective discharge analysis of ecological processes in streams. Water Resour. Res. 41 , W11411 (2005).
4 N. L. Poff , The natural flow regime. BioScience 47 , 769–784 (1997).
5 J. A. Francis, N. Skific, S. J. Vavrus, J. Cohen, Measuring “weather whiplash” events in North America: A new large-scale regime approach. JGR Atmos. 127 , e2022JD036717 (2022).
6 G. Myhre , Frequency of extreme precipitation increases extensively with event rareness under global warming. Sci. Rep. 9 , 16063 (2019).31690736
7 D. B. Booth, C. R. Jackson, Urbanization of aquatic systems: Degradation thresholds, stormwater detection, and the limits of mitigation. J. Am. Water Resour. Assoc. 33 , 1077–1090 (1997).
8 L. E. McPhillips, S. R. Earl, R. L. Hale, N. B. Grimm, Urbanization in arid central Arizona watersheds results in decreased stream flashiness. Water Resour. Res. 55 , 9436–9453 (2019).
9 J. R. Blaszczak, J. M. Delesantro, D. L. Urban, M. W. Doyle, E. S. Bernhardt, Scoured or suffocated: Urban stream ecosystems oscillate between hydrologic and dissolved oxygen extremes. Limnol. Oceanogr. 64 , 877–894 (2019).
10 B. O’Donnell, E. R. Hotchkiss, Resistance and resilience of stream metabolism to high flow disturbances. Biogeosciences 19 , 1111–1134 (2022).
11 A. J. Reisinger , Recovery and resilience of urban stream metabolism following Superstorm Sandy and other floods. Ecosphere 8 , e01776 (2017).
12 B. R. Deemer , Experimental reductions in subdaily flow fluctuations increased gross primary productivity for 425 river kilometers downstream. PNAS Nexus 1 , pgac094 (2022).36741441
13 M. Rode , Sensors in the stream: The high-frequency wave of the present. Environ. Sci. Technol. 50 , 10297–10307 (2016).27570873
14 N. B. Grimm, S. G. Fisher, Stability of periphyton and macroinvertebrates to disturbance by flash floods in a desert stream. J. North Am. Benthol. Soc. 8 , 293–307 (1989).
15 K. Ogle , Quantifying ecological memory in plant and ecosystem processes. Ecol. Lett. 18 , 221–235 (2015).25522778
16 D. Tilman, J. A. Downing, Biodiversity and stability in grasslands. Nature 367 , 363–365 (1994).
17 M. G. Turner, K. H. Braziunas, W. D. Hansen, B. J. Harvey, Short-interval severe fire erodes the resilience of subalpine lodgepole pine forests. Proc. Natl. Acad. Sci. U.S.A. 116 , 11319–11328 (2019).31110003
18 M. L. Kirwan, G. R. Guntenspergen, J. T. Morris, Latitudinal trends in Spartina alterniflora productivity and the response of coastal marshes to global change. Global Change Biol. 15 , 1982–1989 (2009).
19 K. W. Krauss , How mangrove forests adjust to rising sea level. New Phytol. 202 , 19–34 (2014).24251960
20 S. B. Katz, C. Segura, D. R. Warren, The influence of channel bed disturbance on benthic Chlorophyll a: A high resolution perspective. Geomorphology 305 , 141–153 (2018).
21 E. H. Stanley, S. M. Powers, N. R. Lottig, The evolving legacy of disturbance in stream ecology: Concepts, contributions, and coming challenges. J. North Am. Benthol. Soc. 29 , 67–83 (2010).
22 D. C. Reed , Wave disturbance overwhelms top-down and bottom-up control of primary production in California kelp forests. Ecology 92 , 2108–2116 (2011).22164835
23 A. P. Appling , The metabolic regimes of 356 rivers in the United States. Sci. Data 5 , 180292 (2018).30532078
24 J. R. Blaszczak, C. B. Yackulic, R. K. Shriver, R. O. Hall Jr., Models of underlying autotrophic biomass dynamics fit to daily river ecosystem productivity estimates improve understanding of ecosystem disturbance and resilience. Ecol. Lett. 26 , 1510–1522 (2023).37353910
25 U. Uehlinger, Resistance and resilience of ecosystem metabolism in a flood-prone river system. Freshwater Biol. 45 , 319–332 (2000).
26 E. S. Bernhardt , Light and flow regimes regulate the metabolism of rivers. Proc. Natl. Acad. Sci. U.S.A. 119 , e2121976119 (2022).35165178
27 M. Arroita, A. Elosegi, R. O. Hall, Twenty years of daily metabolism show riverine recovery following sewage abatement. Limnol. Oceanogr. 64 , S77–S92 (2019).
28 L. B. Leopold, T. J. Maddock, The Hydraulic Geometry of Stream Channels and Some Physiographic Implications (USGS Numbered Series No. 252) (Professional paper U.S. Government Publishing Office, 1953).
29 P. Savoy, J. W. Harvey, Predicting light regime controls on primary productivity across CONUS river networks. Geophys. Res. Lett. 48 , e2020GL092149 (2021).
30 J. P. Julian, E. H. Stanley, M. W. Doyle, Basin-scale consequences of agricultural land use on benthic light availability and primary production along a sixth-order temperate river. Ecosystems 11 , 1091–1105 (2008).
31 J. C. Finlay, Stream size and human influences on ecosystem production in river networks. Ecosphere 2 , art87 (2011).
32 J. B. Heffernan , Macrosystems ecology: Understanding ecological patterns and processes at continental scales. Front. Ecol. Environ. 12 , 5–14 (2014).
33 W. E. Ricker, Stock and recruitment. J. Fish. Board Canada 11 , 559–623 (1954).
34 N. L. Poff, J. D. Olden, D. M. Merritt, D. M. Pepin, Homogenization of regional river dynamics by dams and global biodiversity implications. Proc. Natl. Acad. Sci. U.S.A. 104 , 5732–5737 (2007).17360379
35 L. Genzoli, R. O. Hall, Shifts in Klamath River metabolism following a reservoir cyanobacterial bloom. Freshwater Sci. 35 , 795–809 (2016).
36 L. Jiménez, A. Freixa, N. Besolí, S. Sabater, Resistance but not recovery is related to the role of specialist taxa in river communities submitted to hydric stress. Sci. Total Environ. 871 , 161952 (2023).36740050
37 N. S. Marzolf, M. Ardón, Ecosystem metabolism in tropical streams and rivers: A review and synthesis. Limnol. Oceanogr. 66 , 1627–1638 (2021).
38 C. E. Schaum , Temperature-driven selection on metabolic traits increases the strength of an algal–grazer interaction in naturally warmed streams. Glob. Change Biol. 24 , 1793–1803 (2018).
39 J. R. Welter , Does N2 fixation amplify the temperature dependence of ecosystem metabolism? Ecology 96 , 603–610 (2015).26236857
40 M. E. Power, M. S. Parker, W. E. Dietrich, Seasonal reassembly of a river food web: Floods, drought, and impacts of fish. Ecol. Monogr. 78 , 263–282 (2008).
41 Y. Vadeboncoeur, M. E. Power, Attached algae: The cryptic base of inverted trophic pyramids in freshwaters. Ann. Rev. Ecol. Evol. Syst. 48 , 255–279 (2017).
42 F. Keck, F. Lepori, Can we predict nutrient limitation in streams and rivers? Freshwater Biol. 57 , 1410–1421 (2012).
43 A. J. Reisinger, J. L. Tank, M. M. Dee, Regional and seasonal variation in nutrient limitation of river biofilms. Freshwater Sci. 35 , 474–489 (2016).
44 E. S. Bernhardt , The metabolic regimes of flowing waters. Limnol. Oceanogr. 63 , S99–S118 (2018).
45 S. G. Fisher, L. J. Gray, N. B. Grimm, D. E. Busch, Temporal succession in a desert stream ecosystem following flash flooding. Ecol. Monogr. 52 , 93–110 (1982).
46 K. Qasem, S. Vitousek, B. O’Connor, T. Hoellein, The effect of floods on ecosystem metabolism in suburban streams. Freshwater Sci. 38 , 412–424 (2019).
47 B. J. Roberts, P. J. Mulholland, W. R. Hill, Multiple scales of temporal variability in ecosystem metabolism rates: Results from 2 years of continuous monitoring in a forested headwater stream. Ecosystems 10 , 588–606 (2007).
48 S. S. Roley, J. L. Tank, N. A. Griffiths, R. O. Hall, R. T. Davis, The influence of floodplain restoration on whole-stream metabolism in an agricultural stream: Insights from a 5-year continuous data set. Freshwater Sci. 33 , 1043–1059 (2014).
49 S. L. Nielsen, S. Enriquez, C. M. Duarte, K. Sand-Jensen, Scaling maximum growth rates across photosynthetic organisms. Funct. Ecol. 10 , 167 (1996).
50 L. Maltchik, F. Pedro, Responses of aquatic macrophytes to disturbance by flash floods in a Brazilian semiarid intermittent stream. Biotropica 33 , 566–572 (2001).
51 K. V. Sarkanen, C. H. Ludwig, Lignins: Occurrence, Formation, Structure and Reactions (Wiley-Interscience, 1971).
52 M. G. Wolman, J. P. Miller, Magnitude and frequency of forces in geomorphic processes. J. Geol. 68 , 54–74 (1960).
53 S. M. Lawlor, “Determination of channel-morphology characteristics, bankfull discharge, and various design-peak discharges in Western Montana” (Scientific Investigations Report 2004-5263, U.S. Geological Survey, 2004).
54 C. I. Mulvihill, B. P. Baldigo, S. J. Miller, D. DeKoskie, J. DuBois, “Bankfull discharge and channel characteristics of streams in New York State” (Scientific Investigations Report 2009-5144, U.S. Geological Survey, 2009).
55 Maryland Department of the Environment, Maryland’s Waterway Construction Guidelines (Maryland Department of the Environment, 2000).
56 L. G. Larsen, J. W. Harvey, Disrupted carbon cycling in restored and unrestored urban streams: Critical timescales and controls. Limnol. Oceanogr. 62 , S160–S182 (2017).
57 B. L. O’Connor, J. W. Harvey, L. E. McPhillips, Thresholds of flow-induced bed disturbances and their effects on stream metabolism in an agricultural river. Water Resour. Res. 48 , W08504 (2012).
58 T. A. Warnaars, M. Hondzo, M. E. Power, Abiotic controls on periphyton accrual and metabolism in streams: Scaling by dimensionless numbers. Water Resour. Res. 43 , W08425 (2007).
59 M. J. Bernot , Inter-regional comparison of land-use effects on stream metabolism: Inter-regional stream metabolism. Freshwater Biol. 55 , 1874–1890 (2010).
60 W. F. Cross , Nutrient enrichment intensifies the effects of warming on metabolic balance of stream ecosystems. Limnol. Oceanogr. Lett. 7 , 332–341 (2022).
61 B. O. L. Demars , Impact of warming on CO2 emissions from streams countered by aquatic photosynthesis. Nat. Geosci. 9 , 758–761 (2016).
62 P. Savoy , Data and code from “Light and flow regimes regulate the metabolism of rivers”. Figshare. 10.6084/M9.FIGSHARE.C.5812160.V4. Deposited 9 March 2023.
63 J. R. Blaszczak , Data from “Extent, patterns, and drivers of hypoxia in the world’s streams and rivers”. ScienceBase. 10.5066/P99X6SIR. Deposited 4 May 2021.
64 L. McKay , NHD Plus Version 2: User Guide. United States Environmental Protection Agency and U. S. Geological Survey. https://www.epa.gov/system/files/documents/2023-04/NHDPlusV2_User_Guide.pdf. Accessed 4 May 2021.
65 R. A. Hill, M. H. Weber, S. G. Leibowitz, A. R. Olsen, D. J. Thornbrugh, The Stream-Catchment (StreamCat) dataset: A database of watershed metrics for the conterminous United States. J. Am. Water Resour. Assoc. 52 , 120–128 (2016).
66 U. S. Geological Survey, USGS water data for the Nation. 10.5066/F7P55KJN. Accessed 3 March 2023.
67 L. A. De Cicco, D. Lorenz, R. M. Hirsch, W. Watkins, M. Johnson, dataRetrieval: R packages for discovering and retrieving water data available from U.S. federal hydrologic web services (Version 2.7.11, 2021).
68 G. M. Hornberger, J. P. Raffensperger, P. L. Wiberg, K. N. Eshleman, Elements of Physical Hydrology (The Johns Hopkins University Press, 1998).
69 J. M. Buffington, D. R. Montgomery, A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers. Water Resour. Res. 33 , 1993–2029 (1997).
70 G. Parker, Surface-based bedload transport relation for gravel rivers. J. Hydraul. Res. 28 , 417–436 (1990).
71 M. Auger-Méthé , A guide to state–space modeling of ecological time series. Ecol. Monogr. 91 , e01470 (2021).
72 B. Rosenbaum, M. Raatz, G. Weithoff, G. F. Fussmann, U. Gaedke, Estimating parameters from multiple time series of population dynamics using Bayesian inference. Front. Ecol. Evol. 6 , 234 (2019).
73 T. Binzer, K. Sand-Jensen, A.-L. Middelboe, Community photosynthesis of aquatic macrophytes. Limnol. Oceanogr. 51 , 2722–2733 (2006).
74 Stan Development Team, Stan Modeling Language Users Guide and Reference Manual (Version 2.26.1, 2022). https://mc-stan.org. Accessed 9 May 2023.
75 R Core Team, R: A Language and Environment for Statistical Computing (R Project for Statistical Computing, 2021).
76 Stan Development Team, RStan: The R interface to Stan (Version 2.26.23, CRAN, 2021). https://CRAN.R-project.org/package=rstan. Accessed 28 September 2023.
77 RStudio Team, RStudio: Integrated Development Environment for R (Version 2023.06.2 + 561, Posit, 2021). https://posit.co/. Accessed 30 August 2023.
78 Advanced Research Computing Center, Beartooth Computing Environment, x86_64 cluster. University of Wyoming, Laramie, WY (2023).
79 H. Wickham , Welcome to the Tidyverse. J. Open Source Softw. 4 , 1686 (2019).
80 K. Müller, here: A Simpler Way to Find Your Files (Version 1.0.1, CRAN, 2020). https://CRAN.R-project.org/package=here. Accessed 5 September 2023.
81 J. Gabry, shinystan: Interactive Visual and Numerical Diagnostics and Posterior Analysis for Bayesian Models (Version 2.6.0, CRAN, 2018). https://CRAN.R-project.org/package=shinystan. Accessed 28 September 2023.
82 S. Xu , Use ggbreak to effectively utilize plotting space to deal with large datasets and outliers. Front. Genet. 12 , 774846 (2021).34795698
83 T. L. Pedersen, patchwork: The Composer of Plots (Version 1.1.3, CRAN, 2020). https://cran.r-project.org/package=patchwork. Accessed 28 September 2023.
84 R. Hilborn, Simplified calculation of optimum spawning stock size from Ricker’s stock recruitment curve. Can. J. Fish. Aquat. Sci. 42 , 1833–1834 (1985).
85 P.-C. Bürkner, brms: An R package for Bayesian multilevel models using Stan. J. Stat. Soft. 80 , 1–28 (2017).
86 A. Vehtari, A. Gelman, J. Gabry, Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27 , 1413–1432 (2017).
87 M. Kay, tidybayes: Tidy data and geoms for Bayesian models (Version 3.0.6, CRAN, 2023). https://CRAN.R-project.org/package=tidybayes. Accessed 28 September 2023.
