Fingerprinting hydrological and biogeochemical drivers of freshwater quality

Understanding the interplay between hydrological flushing and biogeochemical cycling in streams is now possible owing to advances in high‐frequency water quality measurements with in situ sensors. It is often assumed that storm events are periods when biogeochemical processes become suppressed and longitudinal transport of solutes and particulates dominates. However, high‐frequency data show that diel cycles are a common feature of water quality time series and can be preserved during storm events, especially those of low‐magnitude. In this study, we mine a high‐frequency dataset and use two key hydrochemical indices, hysteresis and flushing index to evaluate the diversity of concentration‐discharge relationships in third order agricultural stream. We show that mobilization patterns, inferred from the hysteresis index, change on a seasonal basis, with a predominance of rapid mobilization from surface and near stream sources during winter high‐magnitude storm events and of delayed mobilization from subsurface sources during summer low‐magnitude storm events. Using dynamic harmonic regression, we were able to separate concentration signals during storm events into hydrological flushing (using trend as a proxy) and biogeochemical cycling (using amplitude of a diel cycle as a proxy). We identified three groups of water quality parameters depending on their typical c‐q response: flushing dominated parameters (phosphorus and sediments), mixed flushing and cycling parameters (nitrate nitrogen, specific conductivity and pH) and cycling dominated parameters (dissolved oxygen, redox potential and water temperature). Our results show that despite large storm to storm diversity in hydrochemical responses, storm event magnitude and timing have a critical role in controlling the type of mobilization, flushing and cycling behaviour of each water quality constituent. Hydrochemical indices can be used to fingerprint the effect of hydrological disturbance on freshwater quality and can be useful in determining the impacts of global change on stream ecology.


| INTRODUCTION
High-frequency monitoring of freshwaters with in situ sensors enables collection of hydrochemical data at time scales matching those of hydrological, biological and chemical dynamics observed in streams (Kirchner et al., 2004). This unprecedented technological breakthrough enables us to address fundamental questions about how stream ecosystems function, including process understanding of dominant hydrochemical and stream metabolic regimes (Bernhardt et al., 2018;Burns et al., 2019). Insights regarding the stream hydrochemical regime and the interplay between hydrological and biogeochemical processes are encapsulated in the concentrationdischarge (c-q) relationship and can be interrogated for high and low flows with the aid of high-frequency measurements.
During high flows and under hydrological flushing, the downstream transport of solute and particulate material dominate (Boyer et al., 1997;Burns, 2005), and high-frequency measurements enable the detection of hysteresis responses indicating changes in solute/particulate mobilization and delivery. Several studies have attempted to quantify storm event c-q responses by designing hydrochemical indices to describe the direction of the hysteresis loops and their magnitude; with hysteresis and flushing indices being the most popular (Butturini et al., 2008;Lloyd et al., 2016a). Although several relationships were reported between the c-q indices and potential controlling factors (Aubert et al., 2016;Bieroza & Heathwaite, 2015;Moatar et al., 2017), we lack a systematic understanding of how different hydrological and biogeochemical processes contribute to the observed hydrochemical patterns.
During low flows, when biogeochemical cycling in freshwaters dominates, high-frequency measurements enable detection of diel and sub-daily cycles in hydrochemical time series. Diel cycles are driven by diurnal changes in light, temperature, redox potential, stream metabolism (production and respiration), uptake/release, denitrification/nitrification, evapotranspiration and bioturbation Nimick et al., 2011;Smolders et al., 2017). Owing to the dependency of biogeochemical processes on temperature, diel cycles are most pronounced during spring and summer, with lower amplitudes in winter (Burns et al., 2016;Halliday et al., 2012). The change in amplitude of diel cycles expresses the intensity of biogeochemical processes and is often used to derive an estimation of their rates, for example, for stream metabolism (Bernhardt et al., 2018) and nutrient assimilation (Rode, Halbedel Nee Angelstein, et al., 2016).
The interplay between hydrological and biogeochemical processes at different spatial scales (horizontal from hillslopes to riparian zone and vertical from groundwater to hyporheic zone) is controlled by catchment properties including water travel times and the rates of biogeochemical processes responsible for retention, mobilization and delivery of solutes and particulates to the stream network (Ameli et al., 2017;Ensign & Doyle, 2006;Wollheim et al., 2018). The land-stream interface of the riparian and hyporheic and zones are areas of particularly dynamic and spatially-variable interactions between hydrological and biogeochemical processes (Johnes et al., 2020;Lewandowski et al., 2020). Generally, the presence of an abundant source of solute in the catchment, for example, agricultural land use for nutrients (Banwart, 2011;Heathwaite, 2010) and bedrock for weathering solutes (Beerling et al., 2020), creates an internal loading that reduces the diversity in the stream c-q responses and increases the frequency of chemostatic behaviour with increasing stream order Godsey et al., 2009Godsey et al., , 2019. However, anthropogenic solutes (nutrients) can show highly variable c-q responses from strong dilution to strong concentration, and the inter-catchment variability is often larger than within-catchment temporal variability due to different geological substrates and depths of solute generation (Botter et al., 2020). Much research has focused on headwater catchments because they are thought to be locations in the stream network where hydrological flushing-biogeochemical cycling and chemostaticchemodynamic responses are potentially in balance, according to the River Continuum conceptual framework (Vannote et al., 1980): there is growing evidence from catchment studies to support this hypothesis (Abbott et al., 2018;Bieroza et al., 2018;Creed et al., 2015;Lynch et al., 2019;Wollheim et al., 2018).
Our understanding of temporal drivers of the interplay between hydrological and biogeochemical processes is continually evolving with intense focus on the analysis of changes in seasonal drivers (temperature and storm event magnitude) and inter-storm variation in solute biogeochemical retention and its effect on flushing behaviour (Burns et al., 2019). Seasonality was found to play an important role in controlling the c-q responses of phosphorus and turbidity, with winter chemostatic-clockwise and summer chemodynamic-anticlockwise storm events in a groundwater-fed river system (Bieroza & Heathwaite, 2015). Seasonal c-q patterns were found to differ from event-based ones for phosphorus and nitrate in 219 French catchments covering a wide range of soil, climate and land use characteristics (Minaudo et al., 2019). These differences between long-term and event-based c-q responses likely indicate different controls on solute transport and storage processes on different temporal scales (Knapp et al., 2020). More recently, the notion that biogeochemical cycling is switched off or severely dampened, depending on the magnitude of the storm event has gained prominence, with higher magnitude storm events thought to have a greater damping effect (Bernhardt et al., 2018;Raymond et al., 2016) compared to low-magnitude storm events that can preserve diel cycles (Burns et al., 2016).
Despite this wealth of information on c-q patterns at high and low flows thanks to advances in high-frequency measurements, our understanding of the interplay between hydrological flushing and instream biogeochemical cycling remains limited. The existing studies tend to focus on either flushing during storm events or cycling during low flows. Here, we aim to advance the understanding of the interplay between flushing and cycling during storm events for multiple parameters and hydrological years. We use high-frequency hydrochemical data to evaluate the variability in storm c-q patterns and seek to partition and understand the discrete contributions from hydrological flushing and biogeochemical cycling. For the purpose of this study we assume that these two processes operate separately over the short timescales of individual storm events (here on average 17 hours), and that both contribute to the observed solute and particulate c-q signals.
We base our analysis on the calculation of hydrochemical indices (i.e. the hysteresis and flushing indices) for whole and decomposed (into underlying trend and diel component) concentration time series.
We synthesise our findings into a conceptual framework of solute and particulate hydrochemical behaviour under storm events of varying magnitude for a third order agricultural catchment in a groundwaterfed river system that has been the focus of detailed research on hyporheic zone processes Heppell et al., 2013;Lansdown et al., 2015) and a broader study of the entire river Eden catchment (Ockenden et al., 2016;Reaney et al., 2019).

| Study catchment and hydrochemical data
The River Leith in situ automated laboratory (2009)(2010)(2011)(2012)(2013)(2014) provided measurements of total phosphorus (TP), total reactive phosphorus (TRP), turbidity (TURB), nitrate nitrogen (NO 3 -N), dissolved oxygen (DO), specific conductivity (COND), redox potential (RED), pH and water temperature (TEMP) for unfiltered stream water samples on an hourly basis. TP and TRP measurements were carried out by the wetchemistry analyser (Systea's Micro Mac) with the remaining parameters analysed using Systea's WaterWatch instrument. The details of the experimental setup can be found in (Bieroza et al., 2014(Bieroza et al., , 2018. Continuous 15-min stream stage, flow discharge and meteorological data were obtained from the UK National River Flow Archive (https:// nrfa.ceh.ac.uk/, UK Centre for Ecology and Hydrology).
The River Leith catchment (54 km 2 , third Strahler order, 2. 6 W 54.5 N, rainfall 957 mm, 1999-2014) is dominated by clay loam and silty loam soil textures, is predominantly used for agriculture and has ca. 85% permanent grassland cover. The study site near Cliburn has been the subject of detailed reach scale evaluation of physical hydrology , nitrogen dynamics (e.g., Krause et al., 2009) and water quality with in situ sensing (Bieroza & Heathwaite, 2015, 2016. The River Leith is a sub-catchment to the river Eden basin (2400 km 2 ), where the predominantly chemostatic stream export of nutrients, particularly N, is controlled by the legacy stores within the Penrith Sandstone (Wang & Burke, 2017). Such chemostatic response reflects large mass legacy stores that buffer the variability in concentrations  and integrate, for example, the long-term increases in nitrate concentrations observed in the sandstone aquifer of the Eden catchment. There is potential for legacy N to accumulate in the stream network over time, which may obscure hydrochemical indices with increasing stream order. Previous work suggests that in the third order River Leith, there remains a dynamic equilibrium between chemostatic and chemodynamic responses (Bieroza et al., 2014(Bieroza et al., , 2018Bieroza & Heathwaite, 2015).

| Data processing and analyses
Hourly water quality data were collated, checked for outliers and temperature-corrected. The data are expressed as values at 25 C for DO, COND and RED. The quality of the dataset was assured by cross-checking with routine (weekly) grab samples and also comparison with monitoring data independently collected by the UK Environment Agency (http://environment.data.gov.uk/water-quality).
Flow discharge data (15 min interval) were used to identify discrete storm events in two steps. The first step used identification of events with flow discharge above 0.3 m 3 s −1 , starting at the first time with an increase in flow of 10% in 1 hr and finishing at the first time with a decrease of 1% following the approach of (Dupas et al., 2016).
The second step involved analysis of the flow discharge first derivative to identify rising and falling limbs as periods with flow above the calibrated threshold values of 0.05 and −0.025 ls −2 respectively . The multiple-peak events were identified as single storm events if the flow discharge between individual peaks had not dropped by more than 50% and if the gap between them was less than 3 hrs. For each storm event, a number of characteristics were calculated including duration of the event, duration of the rising and falling limb, maximum flow discharge, flow discharge volume during the event, time from the previous storm event, time to the next storm event and a 5 day antecedent precipitation index (API 5 ; Equation (1), Tables 1 and S1) (Brocca et al., 2008). We classified storm events according to their magnitude using a single Q max threshold of 1.7 m 3 s −1 , resulting in 80 high-magnitude (average Q max = 13.2 m 3 s −1 ) and 74 low-magnitude (average Q max = 0.72 m 3 s −1 ) storm events respectively. We only examined storm events with >70% data completeness for all tested parameters.
To evaluate the effect of each storm event on water quality parameters, each dataset was divided into three periods: (1) before storm event (BF), (2) during storm event (ST) and (3) after storm event (AF), with ST corresponding to the storm event identified as above and BF and AF being actual times from the previous and to the next storm event. In cases where BF and AF were longer than 30 hrs, a fixed time of 30 hrs was used.
Two empirical indices were calculated to describe c-q responses during storm events: hysteresis index (H i ; Equation (2), Table 1 (Lloyd et al., 2016a). In many instances e.g. from 34% for TEMP to 59% for RED (average for all parameters 42%), it was not possible to calculate the exact H i_50 value due to rising and falling limbs covering different flow percentiles. Consequently, we used a mean value of H i across all available flow percentiles (H i_Mean ). The differences between mean H i_50 and H i_Mean were not statistically significant ( Figure S1).
The flushing index describes a relative change in parameter concentrations during the storm event: F i < −0.1 dilution responseconcentrations decrease with the flow, −0.1 < F i < 0.1 neutral response, F i > 0.1 concentration response-concentrations increase with the flow. We observed that F i , since being based on just two points, was sensitive to the presence of noisy data, multiple-peak storm events and diurnal changes in concentrations ( Figure S4) leading to inaccurate representation of flushing behaviour. Therefore, we propose a new calculation of F i_Mean , based on mean concentrations in the period before and during the storm event (Equation (4), Table 1).
The differences between mean F i_Point and F i_Mean were not statistically significant with the exception of TURB for which mean F i_Mean was lower than F i_Point (0.46 vs. 0.55; Figure S5). Based on the combination of H i and F i , nine discrete c-q types were defined to describe the behaviour of solutes and particulates during storm events ( Figure 1).
To investigate the effect of storm events on hydrological flushing and biogeochemical cycling simultaneously, we decomposed the high-frequency time series into the underlying trend as a proxy for hydrological flushing, and amplitude as a proxy for biogeochemical cycling using Dynamic Harmonic Regression (DHR) as described in Equation (5) in Table 1 (Young et al., 1999). The method is typically used to examine the long-term trends, seasonality and short-term fluctuations in time series, including water quality data (Halliday et al., 2012). The method is also suitable for analysing highfrequency water quality time series since it can handle missing data (Young et al., 1999). Here, we conducted the time series decomposition into the underlying trend and short-term components (periodicity corresponding to 24, 12 and 8 hr components) using the CAPTAIN toolbox (Young et al., 2019). The exact periodicity of the cyclic components was determined with the autoregression model; for the diel component it varied from 21 to 27 hrs. The DHR model estimates the parameters recursively (trend's slope, amplitude and phase of the short-term cyclical components) using the Kalman Filter and the Fixed Interval Smoother (Halliday et al., 2012;Young et al., 1999). The parameters are time-varying, which allows examination of their behaviour before, during and after the storm event. Since the storm event time series are trend dominated, there is a risk of spectral overlap between trend and cyclic components and model overparameterisation. To avoid this, visual inspection of spectrum and decomposed trend and cyclical components was carried out and a single value of Noise Variance Ratio (NVR) hyper-parameter was chosen for all cyclical components to avoid model overfitting (Young et al., 1999(Young et al., , 2019. We focus on the trend and diel component as they explained most of the variance (>90%) in the concentration time series and we assume that they represent hydrological flushing and biogeochemical cycling respectively. To quantify flushing and cycling behaviour, we calculated F i for trend (F i_Trend ) and amplitude (F i_Amp ) separately, using Equation (4) in Table 1. Based on this data, we have defined six dominant flushing and cycling responses: (1) concentra- (4) enhancement effect F i_Amp > 0.1, (5) dampening effect F i_Amp < −0.1, and (6) neutral effect (neither enhancement nor dampening) To compare the differences in mean concentrations between the parameters, seasons and storm characteristics (Table S1), a nonparametric analysis of variance was used (Kruskal-Wallis test). All data processing and statistical analyses were carried out in MATLAB version 9.4 (MathWorks, 2018).
T A B L E 1 Definition of the key variables and their equations used in the study

Equations Definition Description
(1) API 5 = kP 2 + k 2 P 3 + k 3 P 4 + k 4 P 5 + P 1 API 5 -a 5 day antecedent precipitation index; k-is a decay parameter here set to 0.95; P i -sum of precipitation for ith day before the storm event (Brocca et al., 2008) (2) (Lloyd et al., 2016b); k-index calculated at 5% flow discharge, between 0 and 1; H i_k -H i calculated for the kth percentile of flow e.g. H i_50 is a H i calculated for the 50th flow percentile; H i_Mean -H i calculated as the mean value of all available values of H i_k , n number of all available H i_k values; Q k -flow calculated for percentile k; C RL_Qk -normalized concentration on the rising limb corresponding to Q k , calculated as a linear interpolation based on normalized flow and concentrations on the rising limb (Q RL and C RL ); C FL_Qk -normalized concentration on the falling limb corresponding to Q k , calculated as a linear interpolation based on normalized flow and concentrations on the falling limb (Q FL and C FL ); Q max -maximum value of Q during storm event; Q min -minimum value of Q; C maxmaximum concentration; C min -minimum concentration. (3) F i_Point -flushing index calculated based on single concentration values, C sconcentration at the peak of the storm event, C b -concentration at the baseflow, beginning of the storm event, C p -the highest of the two values: concentration for the storm event (ST), C BF -mean concentration for the period before the storm event (BF), C p -the highest of the two values: here assumed to be zero due to the short-term duration of individual storm events, Short-cyclical components of period 24, 12, 8 and 6 hr, e t -error component (Young et al., 1999) 3 | RESULTS

| Diversity in c-q responses
The dataset spans over five hydrological years ( Table S1).
For the selected 154 storm events and nine water quality parameters, hysteresis (H i ) and flushing indices (F i_Point and F i_Mean ) were calculated to describe the diversity in storm event c-q responses (Table 2). In general, H i was close to zero for most of the parameters From nine possible c-q types based on the distribution of H i and F i values (Figure 1), the parameters showed a strong preference towards just four types (  Table 2). From a total of 154 storm events, 151 showed a unique combination of c-q types.
There were significant seasonal differences in H i and F i for several parameters ( Figure 2, Table 2 and Figure S2). Anticlockwise summer (H i < 0) and clockwise winter (H i > 0) storm events were typical for TP, TRP, TURB and NO 3 -N. Both DO and RED showed an opposing trend (summer clockwise/neutral and winter anticlockwise response) and COND and pH showed a strong preference towards clockwise events with higher H i for winter storm events. Fewer parameters showed significant seasonal trends in F i , with a stronger concentration response for TURB and stronger dilution for NO 3 -N and COND in winter (Table 2). Summer storm events showed on average ( Figure S3 Several of the parameters showed significant seasonal differences in F i_Trend and F i_Amp (  T A B L E 3 Diversity of storm event c-q responses: Percentage contribution of each c-q response type (as in Figure 1) For NO 3 -N, DO, COND, RED, pH, the opposite was true. We observed a general concentration effect for storm events with a shorter rising limb and lower Q max . Enhanced diel cycles were observed for lower Q max and longer time from previous storm event for nutrients and for higher Q max and shorter time from previous storm event for COND and DO. We further investigated the impact of Q max on the observed patterns by separating patterns in trend and amplitude for high-and low-magnitude storm events ( Figure S6).
The c-q responses during high-( Figure 5) and low-magnitude

| Is the flushing index a correct representation of concentration changes during storm events?
We found that NO 3 -N, DO, COND, RED, pH and TEMP showed stability in the flushing indices with near zero values of F i for different seasons and storm events, indicating that concentrations at the beginning of the storm event (C b ) were similar to concentrations at the peak flow discharge (C s ). As these parameters exhibit strong diel cycling (Table S2) We compared two flushing indices: F i_point calculated from a single concentration at the beginning of the storm event and at the peak flow discharge (Butturini et al., 2008) and F i_Mean based on average concentrations before and during the storm event. We found that there was a good agreement between F i_Point and F i_Mean overall ( Figure S5) and for the high-magnitude storm events but a poorer agreement for the lowmagnitude storm events with noisy data, diel patterns and multiple peaks ( Table 2). An example of this effect can be seen for a lowmagnitude storm event in Figure 6 and Table S3. The readily available source of solutes or particulates (e.g., erosion of river bed sediments) and/or their delivery along relatively short pathways (e.g. from the hyporheic and riparian zones) (Bernhardt et al., 2017).

Negative values of H i (anticlockwise responses) indicate a distant
source, delayed mobilization and/or longer delivery pathways (e.g., from isolated locations in the catchment or along subsurface flow pathways) (Bieroza & Heathwaite, 2015). The F i defines a concentration change between pre-and during the storm event, and therefore reflects differences in concentrations between old, preevent water and new storm event water. Positive values of F i (concentration responses) indicate that storm water has higher concentra- nutrients and sediments c-q responses was driven by the magnitude of storm events were found in other studies (Butturini et al., 2008;Lloyd et al., 2016b).
Fewer statistically significant seasonal patterns were observed in the flushing responses, which may indicate variation in source/delivery patterns from storm to storm as the controlling factor for the flushing behaviour for TP, TRP and TURB (concentration) and dilution response for NO 3 -N. The strongest concentration response was observed for TURB and TP, indicating that a large portion of the phosphorus pool is driven by sediment transport (Lloyd et al., 2016b). Total reactive P showed a mixed dilution-concentration pattern that may indicate a switch between delivery from shallow pathways (as TURB and TP) and dilution of point sources (e.g., sewerage systems) or Prich shallow subsurface pathways. As both point and groundwater sources should become more pronounced during summer when low flows tend to dominate (Charlton et al., 2018), the apparent lack of seasonal differences in TRP flushing suggests that TRP transport is driven by storm to storm retention and flushing mechanisms. For the remaining parameters (NO 3 -N, DO, COND, RED, pH and TEMP), there was a better correlation with the magnitude of the storm event than with seasonality.

| The interplay between hydrological flushing and biogeochemical cycling
To date, most studies have focused on either flushing and storm events (see e.g., (Bieroza & Heathwaite, 2015;Blaen et al., 2017;Boyer et al., 1997) or cycling and diel patterns, see for example (Burns et al., 2016;Ensign & Doyle, 2006). Our approach allowed us to conceptualise the dual flushing-cycling behaviour in relation to the magnitude of the storm events (Table 4)  The predominance of a given flushing pattern is likely to be linked to the major source of the chemical. For example, Knapp et al. (2020) distinguished different c-q behaviour types for solutes of different origin: predominantly concentration behaviour for soil-sourced solutes, predominantly dilution behaviour for groundwater-sourced solutes and a mixed behaviour for solutes derived from different sources. The chemical origin can also play a role in the predominant direction of the c-q relationship on an event-scale. For example, Rose et al. (2018) showed that geogenic solutes (calcium or nitrate) show predominantly a clockwise response, while biologically associated solutes show predominantly an anticlockwise response. These findings corroborate to some extent our results that low-magnitude storm events with more pronounced biogeochemical cycling exhibit mostly anticlockwise responses, while high-magnitude storm events with more pronounced hydrological flushing show mostly clockwise responses. However, this pattern was not supported for two cycling dominated parameters DO and RED, which suggests that c-q relationships for different catchments, storm events and parameters can be highly variable.
Storm events are periods of rapid changes in physical, chemical and biological variables in riverine systems, with the changes in shear stress, nutrient availability, increased turbidity and reduced light penetration being the most critical. This disturbance can either enhance or dampen the biogeochemical processes and the amplitude of diel cycles (Woodward et al., 2016) with a relative balance of flushing and cycling for low-magnitude storms and increased flushing over cycling for high-magnitude storm events (Raymond et al., 2016;Wollheim et al., 2018). High-magnitude storm events often involve stream bed erosion that leads to scouring or burial of the biomass of autotrophic and heterotrophic communities resulting in diminished (Bernhardt et al., 2018) or completely suppressed (Burns et al., 2016)  TEMP) showed consistently reduced amplitudes of diel cycles suggesting that the occurrence of the disturbance rather than its magnitude was a controlling factor. For the low-magnitude summer storm events specifically, their increased occurrence was shown to increase delivery of organic matter and benthic respiration resulting in reduced concentrations and the amplitude of diel cycles (Hutchins et al., 2020). We observed enhanced rather than reduced cycling during high-magnitude storm events for some parameters in our study (COND,. This behaviour may be linked to flushing of hyporheic flow pathways and increased efficiency of aerobic respiration and denitrification that is not observed during low-magnitude storm events (Trauth & Fleckenstein, 2017).

| Advances in understanding hydrochemical processes in-stream systems
The full potential of the new wave of high-frequency measurements and its impact on advancing stream hydrochemistry and ecology (Kirchner et al., 2004) is only just being realised (Benettin & van Breukelen, 2017;Rode, Wade, et al., 2016). Our ability to detect significant patterns in the c-q data and to conceptualise this behaviour improves with the increasing availability of long-term high-quality and high-frequency datasets, where the challenges with water quality sample degradation over time have been overcome. Improved conceptual models of stream functioning are critical to our understanding of the current and future impacts of multiple stressors on stream ecosystems (Birk et al., 2020). Regional (Ockenden et al., 2017) and global (Seneviratne et al., 2012) hydroclimatic simulations indicate that the frequency of extreme hydrological events is rising, including increased storm and drought occurrence. Such events are predicted to increase in frequency, intensity and duration with knock-on consequences for freshwater ecology (Ouellet et al., 2020;Woodward et al., 2016) including under baseflow conditions (Arnell, 2004) and especially for nutrients in groundwater-fed river systems, for example (Trauth & Fleckenstein, 2017 (Birk et al., 2020). The ratio of flushing indices for trend (a proxy for flushing) and amplitude (a proxy for cycling) could be calculated to evaluate a relative importance of hydrological vs. biogeochemical controls. Such an approach could help to test conceptual frameworks of River Network Saturation concept (Wollheim et al., 2018) and pulse-shunt concept (Raymond et al., 2016) over a range of hydrological conditions and stream networks.
We recommend the calculation of flushing indices to provide insight into average concentration characteristics for periods of interest (e.g., before, during and after the storm event) because reliance on point values as suggested by (Butturini et al., 2008)

| Potential limitations
Our approach applied to decompose the water quality signal into underlying trend and amplitude of diel cycle requires further testing.
While the DHR decomposition method is widely accepted (Halliday et al., 2012;Young et al., 1999), it may be less suited to model c-q dynamics for shorter storm events (<24 hrs). We sought to avoid these issues by modelling a longer time series including the time before, during and after individual storm events (>60 hrs). Further, the risk that observed patterns are model artefacts requires evaluation, particularly whether the trend and cyclical components are independent or there is a spectral overlap and thus a correlation between concentration/enhancement and dilution/dampening patterns. Our results indicate, however, that there was a large variation in flushingcycling responses including also concentration/dampening and dilution/enhancement patterns. Another potential limitation is the possibility of hydrological and biogeochemical processes both contributing to trend and amplitude of the cyclical component. To test this, an independent way of decomposing signals into hydrological and biogeochemical components is needed, for example, by the use of hydrochemical models, data mining approaches and in situ estimation of the rates of biogeochemical processes (Aubert et al., 2016). This alternative approach is particularly needed when no apparent diel cycle can be detected in the water quality time series and the signal decomposition used in our study cannot be justified. Limiting nutrients often do not exhibit diel cycles due to rapid uptake ; however, in highly polluted agricultural streams a N and P co-limitation is common along with limitation from other factors e.g. light availability (Jarvie et al., 2018). A lack of diel cycle can also result from different biogeochemical processes producing asynchronous diel cycles (Nimick et al., 2011).
As our results show, both H i and F i indices are potentially highly sensitive to the presence of diel cycles in modelled time series. As the period of diel cycle versus storm event should be random ( Figure S1), one can expect a large variation in both hysteresis and flushing indices. Indeed, the hysteresis index H i showed a large variation in clockwise and anticlockwise responses. However, since similar patterns in H i were observed for flushing, mixed type and cycling dominated parameters in response to common drivers (magnitude of a storm event and seasonality), we conclude that the H i values presented in our study were not biased. Unlike the H i , the flushing index F i showed a very narrow range for the mixed type and cycling dominated parameters as most of the responses were classified as flushing neutral (−0.1 < F i < 0.1). This flushing type is likely related to a chemostatic response to flow for several parameters in our study (Godsey et al., 2019). The concentrations of NO 3 -N, COND, pH and RED can be buffered (Bieroza et al., 2014) by the presence of an abundant source of solute from the Penrith Sandstone aquifer (Wang & Burke, 2017). The parameters controlling stream metabolism (DO and TEMP) are dominated by the biogeochemical cycling and therefore potentially less responsive to hydrological flushing. These results suggest that F i is potentially better suited to representing typically flushing parameters like TP, TRP and TURB with highly episodic response to flow due to anthropogenic pollution Bieroza et al., 2014). Novel hydrochemical indices are needed to capture the flushing dynamics of parameters showing a chemostatic behaviour.

| CONCLUSIONS
As the frequency of extreme hydrological events is rising, presenting major challenges for the future of freshwater ecosystems (Bunn, 2016), understanding the interplay between the hydrological flushing and biogeochemical cycling is of critical importance (Wollheim et al., 2018). Combining hydrochemical fingerprints, as undertaken in our study, adds value to our understanding of the drivers of freshwater quality and their potential to impact stream ecology. As coined by (Rode, Wade, et al., 2016), we are riding the highfrequency data wave where advances in in situ monitoring enable insights that were not previously observable; this opportunity is akin to the growth in understanding observed in, for example, eDNA in freshwater environments (Seymour, 2019). Freshwater ecosystems underpin human survival regarding water security yet the biodiversity on which their functioning depends remains threatened worldwide despite efforts during the decade of 'Water for Life' (Dudgeon et al., 2006); we are now looking at emergency recovery plans (Tickner et al., 2020). Our capacity to describe environmental impacts (problem identification) and understanding mechanisms, as this paper shows, remains strong. The challenge is applying this understanding towards preventing impacts through predictive science and appropriate regulatory frameworks to deliver sustainable outcomes. To address this research and science-application gap, we suggest piloting hydrochemical indices as metrics that capture key parts of the biogeo-

SUPPORTING INFORMATION
Supporting tables include information on storm events (Table S1), diel cycles (Table S2), flushing and cycling for two storm events of contrasting magnitude (Table S3) and analysis of variance in flushing and cycling parameters (Table S4). Supporting figures include information on difference between H i_50 and H i_Mean ( Figure S1), seasonal differences in H i and F i ( Figure S2), storm characteristics ( Figure S3), the effect of diel cycle on H i and F i ( Figure S4), differences between F i_Point and F i_Mean ( Figure S5) and differences in trend and amplitude ( Figure S6).

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.