Resistance and resilience of stream metabolism to high flow disturbances

Streams are ecosystems organized by disturbance. One of the most frequent and variable disturbances in running waters is elevated flow. Yet, we still have few estimates of how ecosystem processes, such as stream metabolism (gross primary production and ecosystem respiration; GPP and ER), respond to high flow events. Furthermore, we lack a predictive framework for understanding controls on within-site metabolic responses to flow disturbances. Using five years of high-frequency dissolved oxygen data from an urbanand agriculturally-influenced stream, we estimated daily GPP and ER and analyzed 5 metabolic changes across 15 isolated high flow events. Metabolism was variable from day to day, even during lower flows. Thus, we calculated metabolic resistance as the magnitude of departure from the dynamic equilibrium during antecedent lower flows and quantified resilience from the days until GPP and ER returned to the range of antecedent dynamic equilibrium. We evaluated correlations between metabolic resistance and resilience with characteristics of each high flow event, antecedent conditions, and time since last flow disturbance. ER was more resistant and resilient than GPP. GPP was typically suppressed 10 following flow disturbances, regardless of disturbance intensity. In contrast, the ER magnitude of departure increased with disturbance intensity. Additionally, GPP was less resilient and took longer to recover (0 to >9 days, mean = 2.2) than ER (0 to 2 days, mean = 0.6). Given the flashy nature of streams draining human-altered landscapes and the variable consequences of flow for GPP and ER, testing how ecosystem processes respond to flow disturbances is essential to an integrative understanding of ecosystem function. 15

minute intervals by an in situ YSI 6920V2 sonde (Hession et al., 2020;O'Donnell and Hotchkiss, 2019). Because a freeze event impaired DO measurements from the YSI sonde, we gap-filled with calibration-checked and comparable data from an adjacent PME MiniDOT from 2017-09-01 to 2018-04-14 ( Figure A2; O'Donnell and Hotchkiss 2019). We obtained light data and barometric pressure measurements from a nearby weather station. A Campbell Scientific CS451 pressure transducer recorded stage measurements every 10 minutes. Velocity and width measurements were taken over multiple years to create site-specific relationships between stage, velocity, wetted width, and discharge. A stage-discharge relationship was created in 95 2013 and updated in 2018. Sensors were calibrated every 2-4 weeks.
To remove lower-quality sensor data due to sensor error or periods of low flow, we used data cleaning and quality checks as in O'Donnell and Hotchkiss (2019). Briefly, we excluded values below the 1% and above the 99% quantile for physicochemical parameters that were heavily skewed (i.e., turbidity and conductivity). We removed physicochemical values we knew to be unreasonable (e.g., turbidity was cut off at zero). We calculated daily medians of physicochemical parameters for all days that 100 had at least 80% of measurements over the course of the day after confirming the 80% cutoff as one that would not bias daily medians from dates without gaps in sensor measurements. Data from lower flow periods when individual sensors may have been out of water (Hession et al., 2020) were excluded when values were out of range of grab sample calibration checks.

Ecosystem metabolism
We estimated GPP, ER, and K (air-water gas exchange) from diel O 2 (DO), light, and temperature sensor data using the same 105 inverse modeling approach and data as O'Donnell and Hotchkiss (2019). We selected the streamMetabolizer R package for our analyses (Appling et al., 2018a), which uses Bayesian parameter estimation and a hierarchical state space modeling framework to generate daily estimates of GPP, ER, and K that create the best fit between modeled and observed DO data (Appling et al. 2018b; Equation 1; Table 1). 110 We used most of the default model specs for streamMetabolizer. Model convergence was visualized via traceplot in the rstan package (Stan Development Team, 2019) to identify the proper number of burn-in steps (500); we saved 2000 Markov chain Monte Carlo (mcmc) steps from four chains after burn-in. We modeled GPP, ER, and K with both observation error and process error. We calculated mean daily standard error derived from the mcmc-derived distributions of GPP and ER. Additionally, to decrease the chances of equifinality between GPP, ER, and K estimates (Appling et al., 2018b), we constrained day-to-day mcmc steps (2000) specified for our model. Additionally, to avoid using biased estimates of metabolism, we removed K values below the 1% (< 3.38 d -1 ) and above the 99% (> 27.21 d -1 ) quantile of model estimates. 246 days of metabolism estimates were ultimately removed due to these QA/QC criteria, resulting in 1375 days (of 1621 total from 2013-01-08 to 2018-04-14) quality-checked GPP and ER for further analyses.

130
To identify flow events for our analyses of metabolic resistance and resilience, we calculated the percent change in cumulative daily discharge (Q) relative to the day prior (Equation 2).

%∆Q =
We selected isolated flow events that had a greater than 50% Q change relative to the antecedent cumulative daily Q. We defined isolation as a period of three days before and three days after a high flow event where no other flow events exceeding 10% Q 135 change occurred. In total, there were 15 isolated flow events across all 5 years that met our criteria for isolated flow events and had quality-checked metabolism estimates ( Figure 2). Hydrograph and metabolism time series for each isolated flow event are available in Appendix Figures 4 -18.

Characterizing metabolic resistance and resilience
To acknowledge the pulsing, day-to-day variability of GPP and ER, we used metabolism estimates from three days prior to 140 define a range of antecedent metabolism for each isolated flow event. We quantified metabolic responses to flow disturbances by comparing the pre-event metabolic range with event and post-event metabolism rates. To assess resistance, we estimated metabolic magnitude of departure (M) during events to quantify the resistance of GPP and ER to disturbance. We calculated M per isolated flow event by comparing the difference between GPP and ER to the nearest value of the antecedent range (Equation where X event is either GPP or ER (g O 2 m -2 d -1 ) on the day of the isolated flow event. X prior is the maximum or minimum value of GPP or ER from the antecedent range, depending on if the isolated flow event resulted in a stimulated (increased) or repressed (reduced) metabolic response. For instance, if GPP was repressed during a flow event, M was calculated as the difference between GPP for the isolated flow event and the minimum value from the antecedent 3-day range ( Figure 3). If the 150 estimate of GPP or ER on the event day did not fall outside of the antecedent range, the magnitude of departure was zero, thus indicating high resistance. A negative magnitude of departure (M) represents a suppression, and a positive M a stimulation, of GPP or ER relative to the antecedent equilibrium. To quantify the resilience of GPP and ER, we estimated recovery intervals (RI) by counting the number of days until metabolic rates returned to within the range of pre-event values, signifying a return to antecedent dynamic equilibrium (Fig-155 ure 3). If metabolism from the isolated flow event did not fall outside of the antecedent range and M is zero, there was no RI (metabolism cannot recover if it never shifts outside the range of normality). To ensure additional flow events did not obscure the recovery interval of GPP or ER, we stopped counting RI the day before the next event (i.e., if another flow event happened four days later, we stopped counting RI at 3 days), and have noted this in our results as days+ and used different symbols in data figures. To test for statistically significant differences between ER and GPP recovery intervals (RI ER and RI GPP ) and ER 160 and GPP magnitude of departure ( M ER and M GPP ), we ran Welch's t-tests in R (R Core Team, 2018).

Testing controls on metabolic resistance and resilience
Quantifying how different antecedent conditions induce variable responses from GPP and ER is critical to furthering our understanding of stream ecosystem responses to flow disturbances. We assessed three categories of potential predictors of metabolic resistance and resilience: antecedent conditions, characteristics of the isolated flow event, and characteristics of the We visually identified the most recent flow event (% change in cumulative daily discharge > 50) prior to each isolated flow event. We ran bivariate correlation analyses to quantify the strength and directions of linear relationships between predictor variables and metabolic resistance and resilience using the R cor.test function (R Core Team, 2018). All modeling and analyses were conducted in R (R Core Team, 2018).

Flow and metabolism
Stroubles Creek is a hydrologically flashy stream, with frequent high flow events ( Figure 2). Cumulative daily discharge for days with quality-checked metabolism estimates ranged from 66 to 114,408 m 3 d -1 , with a median of 6,230 m 3 d -1 . The 15 isolated flow events selected for analyses were within the mid-high range of all cumulative daily discharge values, and were Both GPP and ER typically recovered from flow-related stimulation or reduction in less than three days (Table 3). There 200 were many isolated flow events where GPP took multiple days to recover but ER never departed from the antecedent dynamic equilibrium (i.e., M = 0; Figure 5). When M GPP and M ER were both greater than zero, ER almost always recovered faster than GPP. RI GPP ranged from 1-9+ d, with an average of 2.5 d (Table 3). RI ER ranged from 1-2 d, with an average of 1.1 d. There was only one isolated flow event where GPP recovered before ER. While ER always recovered before another flow event occurred, GPP did not recover before another flow event for 4 of 15 analyzed events. Excluding the 4 of 15 events when GPP did not 205 recover, recovery intervals for GPP and ER were not correlated ( Figure 5) or significantly different (p = 0.12, α = 0.05).

Controls on metabolic resistance and resilience after a flow disturbance
Although GPP and ER are linked processes, the variables that were moderate or strong predictors of resistance or resilience (r > 0.5) differed between ER and GPP (Table 4; Figure 6). There were no predictors with moderate or stronger relationships for both M GPP and RI GPP . Because the median RI ER was zero, bivariate correlations could not be used to determine potential 210 predictors of ER resilience. Peak discharge of the flow event was negatively correlated with M ER (r = -0.59). The magnitude of each disturbance, characterized by the % change in cumulative daily discharge, was negatively correlated with M GPP (r = -0.46, p = 0.08) and M ER (r = 0.66, p = 0.01, α = 0.05) (Figure 7), but was not significantly correlated with RI GPP (p = 0.54; α = 0.05). Overall, there were multiple environmental controls on metabolic resistance or resilience that were strongly correlated with either GPP or ER, but no significant drivers of both GPP and ER resistance and resilience.

Metabolic resistance and resilience
GPP and ER responded differently to isolated flow events in a heterotrophic stream draining a heterogeneous urban-agricultural landscape. Notably, ER was more resistant than GPP (Figure 1). Of the fifteen high flow events analyzed here, three events stimulated ER and there were two instances of minor GPP stimulation ( Table 3). The flow disturbance of heterotrophic ac-220 tivity was likely balanced by increased inputs of organic carbon from terrestrial sources that stimulated respiration (Roberts et al., 2007;Demars, 2019). We found similar differences in how GPP and ER change with discharge across all days with metabolism estimates from Stroubles Creek (instead of event-specific analyses presented in this paper); GPP decreased but ER was relatively constant on days with higher than median flow (O'Donnell and Hotchkiss, 2019). How often and when GPP and ER respond similarly to flow disturbances may differ among ecosystems as a function of their metabolic balance (GPP:ER), 225 nutrient limitation status, and history of flow disturbance. Ultimately, flow-induced changes disproportionately disturbed GPP relative to ER.
ER was more resilient than GPP. Differences in ER and GPP resilience were likely a result of flow-induced changes to physicochemical parameters (e.g., turbidity) that can also enhance the physical disturbance of flow on GPP (O'Donnell and Hotchkiss, 2019). For instance, sustained periods of high turbidity following a flow disturbance can prolong the recovery of 230 GPP by inhibiting light attenuation (Blaszczak et al., 2019). In contrast, higher resilience of ER is likely a function of greater resistance of ER to disturbances (i.e., smaller M; Table 3) as well as flow-induced ER stimulation. The correlation of M GPP and M ER , but a lack of correlation between RI GPP and RI ER (Figure 5), suggests GPP and ER were temporarily decoupled while recovering, despite similar initial responses of GPP and ER to flow disturbances.
The dynamic nature of stream metabolism, even during low flow periods, must guide how we quantify metabolic responses 235 to disturbance. When we quantified resistance as a deviation from an antecedent average (e.g., Reisinger et al. 2017;Roley et al. 2014), few flow events yielded 100% resistance and every metabolism estimate was classified as a reduction or stimulation.
Without acknowledging the "pulsing steady state" of stream metabolism, we may overestimate disturbances in ecosystem function. In assessing metabolic responses and recovery from smaller flow events relative to the dynamic equilibrium of metabolism at baseflow, we found some of the shortest metabolic recovery intervals recorded in the literature (Figure 8; Table A1). Incorpo-240 rating the pulsing equilibrium of metabolism and standardizing calculations of metabolic recovery dynamics will enable more robust, cross-site comparisons of complex ecosystem response to changes in flow.

Controls on metabolic resistance and resilience after a flow disturbance
While GPP responded similarly to flows regardless of magnitude, ER was more resistant to smaller magnitude isolated flow events. Our hypothesis that isolated flow events of greater magnitudes (i.e., larger % change in cumulative daily discharge) of autotrophs and heterotrophs (Uehlinger, 2000(Uehlinger, , 2006. Primary producers reside in exposed areas on the streambed to access light required for photosynthesis, and are thus more vulnerable to scour than heterotrophic biofilms tucked within, and protected 250 by, substrates in the streambed, sediments, and hyporheic zone (Uehlinger, 2000). At some threshold of higher flows that disturb more protected areas within and below streambeds, we expect ER will decline as flow-induced stress exceeds flow-induced carbon and nutrient subsidies. Future analyses that include event duration may also provide new insights into flow-metabolism dynamics: Do sustained, higher flows change GPP and ER in the same way as a more instantaneous, intense flow event? As is common of long-term characterizations of metabolism in streams, many high flow days had metabolism model outputs that did Environmental conditions on the day of isolated flow events that promote biomass growth, such as high light and temperature, were not significant predictors of ER or GPP recovery intervals. Metabolic recovery trajectories often increase with temperature and PAR (Uehlinger and Naegeli, 1998;Uehlinger, 2000), and consequently may also display a seasonal component, with faster 270 recoveries in spring and slower recoveries in winter (Uehlinger, 2000(Uehlinger, , 2006. We found a negative relationship between RI GPP and both temperature and season ( Figure A20), suggesting that under colder conditions and in winter months, GPP takes longer to recover. While we did not find strong predictors of RI ER among the environmental variables in our dataset, changes in the source, magnitude, and biological lability of organic matter inputs may alter RI ER (Roberts et al., 2007). How do changes in nutrients and organic matter govern the resistance and resilience of GPP and ER? Combining high-frequency nutrient and 275 organic matter quality measurements with metabolic resistance and resilience estimates will offer an improved understanding of how changing resources mediate metabolic responses to flow changes.
We did not find evidence to support a single subsidy-stress response of metabolic changes during and after high flow events ( Figure 1). ER was more resistant than GPP to most flow disturbances (H1). At small-intermediate sized flow disturbances, the response of metabolism was variable (H2,H3), with the greatest range of metabolic stimulation or reduction (i.e., subsidy 280 or stress) observed at smaller flow changes (Figure 7). ER and GPP did not increase or decrease relative to their antecedent "pulsing equilibrium" during several high flow events (H0). Light, temperature, or turbidity may control the variability in metabolic response to smaller flow disturbances. With increasing intensity of flow disturbance, stress and replacement may indeed scale with intensity (H3).

285
Metabolic regimes are punctuated by high flow events that create frequent pulses of stimulated or repressed GPP or ER (e.g., Uehlinger 2006;Beaulieu et al. 2013;Bernhardt et al. 2018). As such, changes in flow play an influential role in the trends and variability in daily metabolism. While geomorphology and disturbance regimes may control metabolic resistance across sites (Uehlinger, 2000;Blaszczak et al., 2019), within-site variability of M and RI may be controlled by the characteristics of each flow event. Differences between ER and GPP response and recovery to flow disturbances at our study site were controlled by 290 higher resistance and resilience of ER relative to GPP. Within this study, our prediction that ER would be more resistant than GPP to flow disturbances was supported, as ER frequently did not even deviate from the antecedent equilibrium. However, ER had less resistance to events of greater magnitude: M ER was negatively correlated with the % change in discharge of flow event, whereas M GPP was not, suggesting that GPP responded similarly to changes in flow regardless of flow magnitude. Metabolic response to small and intermediate flow disturbances was variable: GPP and ER were both stimulated and suppressed. We 295 suggest there may be a resistance threshold to flow disturbances, where controls other than flow magnitude (e.g., season, light, turbidity) might regulate metabolic responses to lower flow changes. Using segmented process-discharge relationships to quantify a resistance threshold of processes to flow disturbances (O'Donnell and Hotchkiss, 2019) may support a more predictive understanding of metabolic response to flow disturbances.
Ultimately, we are entering an era of metabolic data opportunity (e.g., Bernhardt et al. 2018). As time series of metabolism 300 lengthen and modeling tools improve, we envision exciting opportunities to better assess the consequences of isolated flow events as well as the impacts of multiple, sequential high flow disturbances that did not meet our criteria for analyzing isolated flow events in this paper. While the short time periods between high flow events in many streams and rivers make isolating and quantifying functional resistance and resilience an ongoing challenge, including dynamic flow in our assessment of metabolic regimes is a critical next step toward a more holistic understanding of frequently disturbed ecosystems.

305
Code and data availability. All data and results are included in the Appendices of this paper. Supplemental data files and metadata (not  Antecedent median photosynthetic active radiation (µmol m -2 s -1 ) -0.31 0.12 n/a 0.04 Antecedent median discharge (m -3 s -1 ) 0.08 0.34 n/a 0.40 Antecedent median water temperature ( • C) -0.32 -0.14 n/a -0.29 Antecedent median turbidity (nephelometric turbidity unit, NTU) -0.23 0.13 n/a 0.00 Table A1. Literature review of published reduction and recovery intervals (RI) of stream gross primary production (GPP) and ecosystem respiration (ER) after high flow events. If not enough information was given to calculate reduction or RI, we listed as "n/a". *=approximated days of recovery from figure in publication. **=approximation given in publication.