Drastic reduction of nutrient loading to a reservoir alters its resistance to impacts of extreme climatic events

By perturbing ecosystems, extreme climatic events (ECEs) can impair ecosystems’ resistance and resilience to other pressures, leading to cascading effects on the continued provision of their ecosystem services. In aquatic ecology, most of the studies linking impacts of perturbations on ecosystems are based on controlled experiments and modeling, rather than real-world data. Using a 55 year dataset of hydrometeorological and reservoir water quality variables from the Ter catchment in Spain, we fill this gap by applying non-linear dynamics and extreme value theory concepts to test whether trophic state modulates reservoir ecosystem’s response to ECEs. We show that both Granger causality between hydrometeorological and water quality variables and effects of ECEs on reservoir water quality diminish after drastic reduction in nutrient loading, supporting our hypothesis that the ecosystem’s trophic state modulates its resistance to ECEs. Thus, by safeguarding reservoirs from nutrient pollution, water resources managers can ameliorate impacts of ECEs on ecosystem health.


Introduction
Eutrophication, i.e. the over enrichment of aquatic ecosystems with nutrients (Carpenter 2005), is one of the major challenges affecting the provision of aquatic ecosystem services. In its fifth global environment outlook report, the United Nations Environment Programme estimates that over 40% of the world's water bodies are experiencing moderate to heavy eutrophication (Xia et al 2016). Eutrophication increases drinking water treatment costs through increased organic matter from algal blooms (Foley et al 2012) or high iron and manganese concentrations (Munger et al 2019) promoted by anoxia (Carpenter 2005) and threatens water-based recreational activities from the proliferation of toxin producing cyanobacteria blooms (Padedda et al 2017).
Climate change might worsen eutrophication through increasing water temperature, intense or reduced precipitation, changes in wind speed and direction and solar radiation (Xia et al 2016). As the climate is changing, so too are the occurrence, intensity, magnitude and variability of extreme climatic events (ECEs) such as heatwaves, droughts and heavy precipitation (IPCC 2014). In ecology, ECEs are episodes where unusual values of climate variables result in responses of an ecosystem that are outside the range of normal variability (Smith 2013). By perturbing lake ecosystems, ECEs may impede the provision of aquatic ecosystem services.
The resistance and resilience of lakes to impacts of ECEs may depend more on preceding lake conditions than on characteristics of the extreme events themselves (Thayne et al 2021). Resistance is how (in)susceptible the system is to perturbation (Mitra et al 2015, Nimmo et al 2015, whereas ecological resilience is the extent to which a system recovers (to its original or alternative state) after a perturbation (Perfecto et al 2019, Thayne et al 2021, or the amount of disturbance a system can withstand without flipping into another state (Folke et al 2004).
The degree to which ECEs impact eutrophication related processes in water bodies varies with flow rate, morphology, geographical location and season, making the generalizations of responses impractical (Delpla et al 2009, Mosley 2015. For instance, while high air temperatures have been reported to increase stratification, leading to eutrophication in lakes from internal nutrient loading (Collins et al 2019), elsewhere, stratification is reported to have reduced nutrients in surface water layers (Xia et al 2016). Droughts have been reported to increase eutrophication and algal blooms in mesotrophiceutrophic lakes (Lisi and Hein 2019), due to reduced lake volumes, which concentrates nutrients and promotes cyanobacteria blooms (Wright et al 2014); and increased frequency of water column mixing events (Soares et al 2019), which promotes internal nutrient loading, hence fuelling cyanobacterial blooms. Yet, in the Mediterranean region, lower water levels are reported to have prevented cyanobacteria blooms by promoting the abundance of macrophyte that outcompeted the bloom causing cyanobacteria (Bakker and Hilt 2016). While intense rainfall has been reported to decrease algal blooms due to the reduction of light intensity through high turbidity and lake flushing (Xia et al 2016), in the floodplain lakes of the lower Dutch Rhine river, high rainfall intensity caused cyanobacterial blooms due to an influx of high nutrient river water through floods (Bakker and Hilt 2016). Also, floods are reported to have caused increased phytoplankton production in Gollinsee Lake due to anoxia-induced phosphorus release from the sediments (Brothers et al 2014). Furthermore, for shallow lakes, rising water levels inhibit light transmittance to lake bed macrophytes, leading to loss of their competitive advantage over phytoplankton and the subsequent dominance of cyanobacterial species (Bakker and Hilt 2016).
Such varied responses call for identifying modulators of lake responses to climate extremes, and we posit that lakes' trophic state might be one of the key factors.
The transient nature of most ECEs prompts researchers to rely on information, automatically collected at high frequency in lakes and rivers, to study their impacts on water quality (Marcé et al 2016). However, such a reliance on high frequency information is, to some extent, due to lack of appropriate methods to address these questions using low frequency (e.g. monthly) data, thereby side-lining the much more abundant long-term water quality records gathered from many aquatic ecosystems across the world. In addition, the majority of studies linking impacts of ECEs on ecosystems arise from highly controlled experiments (Veraart et al 2011) and modelling (Dakos et al 2012), hence are far removed from real-world data (Veraart et al 2011). In this study, an innovative time-series approach is applied to long term monthly monitoring data, from a real-world system, to assess the role of the ecosystem's trophic state in modulating the resistance of water quality to impacts of ECEs. We took advantage of data from the long-term water quality monitoring program at Sau Reservoir (Spain) to unravel the behaviour of water quality under ECEs. During the last 55 years, Sau switched between two contrasting trophic states, following a drastic reduction in nutrient loading. We tested the hypothesis that the change in the trophic state of the reservoir significantly increased the resistance of the system's water quality to impacts of ECEs.

Study site
Sau Reservoir (41 • 58 ′ N; 2 • 22 ′ E, maximum depth 60 m, area 7.6 km 2 , mean hydraulic retention time of 90 days) is located on the Ter River, in North East Spain. Its catchment area is 1790 km 2 , with forestry (78%), arable farming (16%) and urbanization as the main land uses (López et al 2011). Sau Reservoir was built in 1963 for hydropower and water supply to 4 million people in the Barcelona Metropolitan area (Ordóñez et al 2010). Sau Reservoir lies downstream from a densely populated area where, between the years of 1991 and 1995, 16 wastewater treatment plants (WWTPs) were built, designed to treat waste water for 530 563 inhabitant-equivalents, representing 95% of the population upstream from the reservoir. The WWTPs were upgraded to secondary treatment during the early 2000s (Marcé et al 2006). The eutrophication process of Sau Reservoir is well described in the literature (Armengol et al 1986, Vidal and Om 1993, Armengol et al 1999, Marcé et al 2004, 2008b, 2010, showing an increase in nutrients loads from the 1960s to the early 1990s, accompanied by an increase in pigment concentrations and anoxia extent. From the late 1990s onwards, nutrient loads and anoxia have decreased substantially as a result of the WWTPs built upstream from the reservoir, which entailed a major reduction of organic matter and nutrient loads to the reservoir, with a sharp decrease of ammonium to nitrate molar ratio, from ∼2.5 to 0.1 (Marcé et al 2008b). These changes improved the trophic status of the system, which previously experienced recurrent episodes of algal blooms and extended deep water anoxia (Ordóñez et al 2010).

Data
This study utilized (a) daily 20 km grid interpolated rainfall (Kg m −2 s −1 ) and air temperature ( • C) of the Escenarios-PNACC dataset (Herrera et al 2012(Herrera et al , 2016) from 1950 to 2015; (b) daily observations of hydrological time series of inflows (m 3 s −1 ) and reservoir water level (m.a.s.l) and monthly reservoir water quality records of temperature (T w , • C), dissolved oxygen (mg DO/l), nitrate (mg NO 3 − /l), ammonium (mg NH 4 + /l) and total phosphorus (mg TP/l) from 1964 to 2019, sampled from several depths at the deepest part of the reservoir. The data was supplied by the Ens d'Abastament d'Aigua Ter-Llobregat. Further details on the water quality data are explained elsewhere in Marcé et al (2008b), López et al (2011) and Šimek et al (2011).

Heatwave indices
Computation of heatwave indices applied the Heatwave Magnitude Index Daily (HWMID) framework (Russo et al 2015) by transforming the time series of daily maximum air temperature (Tmax), from Escenarios-PNACC dataset, into an index. In this framework, a heatwave is defined as a period of 3 consecutive days in which the maximum temperature is above the 90th percentile of the daily maxima in the reference period of 1981-2010, centred on a 31 day window (text S2). HWMID calculations were implemented in the 'extRemes 2.0-0' package (Gilleland and Katz 2016).

Drought/wetness indices
Drought and wetness indices were derived from precipitation, streamflow and reservoir water level time series, using the non-parametric Standardized Drought Analysis Toolbox (Hao and AghaKouchak 2014, Farahmand and AghaKouchak 2015) (text S1), resulting into the non-parametric standardized streamflow index (nSSFI), the non-parametric standardized water level index (nSLI), and the nonparametric standardized precipitation index (nSPI).

Epilimnion and hypolimnion water quality time series
To simplify the analyses, multi-depth water quality time series of T w , DO, NO 3 − , NH 4 + , and TP were aggregated into an epilimnion and a hypolimnion time series. To achieve the aggregation, we supplied the LakeAnalyzer code (Read et al 2011) with a time series of water temperature and reservoir depth profiles, from which the water temperature data was measured, for each day of the record. The code used such data to generate the top (metaT) and bottom (metaB) depths (m). The metaT is essentially the end of the epilimnion, whereas the metaB is the beginning of the hypolimnion. Thus, for water temperature and other reservoir water quality variables recorded at the same depth profiles supplied to the numerical code, the average value, on a particular day, was obtained by calculating the mean of all values at depths equal to and above metaT (for the epilimnion) and all values at depths equal to and below the metaB (for the hypolimnion). The output of this analysis was a single value of water quality variable for the epilimnion and hypolimnion for a particular day. Missing data in the monthly time series, for all water quality variables, were filled by interpolating small gaps of up to 3 months.

Statistical analyses 2.4.1. Granger-causality between hydrometeorological and reservoir water quality variables
This analysis aimed at establishing evidence of causality between the hydrometeorological and epilimnetic & hypolimnetic water quality timeseries, before and after WWTPs. The nonparametric Granger causality in quantiles method (Balcilar et al 2017), a non-linear variant of the Granger causality framework (Granger 2008), was applied to derive causal evidence. In the Granger framework, a causal model between two linear stationary time series X t and Y t is expressed as: where η c is an uncorrelated white-noise series in the complete model. In the model, X t is assumed to cause Y t only if it contains information in previous terms that improves the prediction of Y t and that information is not present in any other predictor (Granger 1988 (1)) brings predictability beyond the one provided by the autocorrelation of the dependent variable at lag 1 and beyond, in the restricted model: where η r is the uncorrelated white noise series of the restricted model. Jeong et al (2012) extended the Granger framework to analysing causality in different quantiles of a time series, because other studies showed evidence of Granger-causality in tails of a variable distribution even when the mean or median did not show any evidence. They derived a 'distance measure' test statistic (J) for testing Granger causality in quantiles of a dependent variable (Fan and Li 1999) based on the comparison of the quantiles in the prediction of Y t with and without considering previous values of X t .
In this study, the Jeong et al (2012) methodology was used to test for evidence of causality between hydrometeorological variables as predictors and epilimnetic and hypolimnetic reservoir water quality time series as response variables, using one predictor and one response variable at a time, before and after WWTPs. Granger causality was calculated at 15 discrete equidistant quantiles between 0 and 1. Significance was set at p ∼ 0.001 for each analysis, running 1000 instances of the calculations with a random predictor, and comparing this outcome with the actual distribution of J values across quantiles ( figure S2(a)). The J statistic, between the predictor and the response variable, is evaluated at 15 quantiles along the range of 0-1, and compared to the area occupied by results coming from 1000 realizations of the same calculation that uses a random predictor instead. We consider occurence of Granger causality, in a given quantile, when the J statistic calculated with the observed predictor is higher than the maximum of the 1000 values obtained with the random predictor, in the same quantile ( figure S2(a)).
In order to synthesize all the information generated by this analysis (1200 quantiles tested across layers, periods, predictors, and response variables), the J statistic was aggregated over three quantile groups (low quantiles: <0.25; central quantiles: 0.25-0.75; high quantiles: >0.75). This was done by calculating the area in the J statistic plot (figure S2(b)), contained between the J statistic line of the predictor and the response variable and the maximum values of the 1000 realizations with the random predictor. Calculated over the three quantile ranges mentioned above, this procedure results in three areas, one for each quantile range (colored areas in figure S2(b)). These three areas were then normalized by dividing each area by the area under the line defined by the maximum values of the 1000 realizations obtained with the random predictors (grey area in figure S2(b)).

Impacts of ECEs on water quality before and after WWTPs
Impacts of ECEs on water quality were assessed by comparing median values of water quality variables in extreme and non-extreme conditions. The heatwave and drought/wetness indices (see sections 2.3.1 and 2.3.2) were used for the identification of ECEs. For the hydrometeorological events, first, nSSFI, nSLI, and nSPI index values considered beyond normal conditions (below or above zero), were clustered using the peaks-over-threshold (POT) approach (Ghil et al 2011), to remove autocorrelation. Then, the distribution of water quality values associated with the ECEs was split twice; first, into periods of before and after WWTPs and, second, corresponding to events with indices beyond the threshold that could be considered extreme (<−1.3 & >+1.3) and the rest that could be considered beyond normal but not extreme (>−1.3 <0 & >0 < 1.3). Thereafter, the value of the water quality variables in the epilimnion and hypolimnion, for each event, was calculated as the median of all data corresponding to that type of event. POT analysis was performed using the POT package, version 1.1-7 (Ribatet 2011). The medians of water quality variables, in extreme and non-extreme events, for each period, were compared using the Welch's test (Zheng et al 2013) on ranked values, which allows a robust comparison of non-normal distributions with different sample sizes and variances. The results of the Welch's test were reported as p-values, that were supplemented by the Hedges's g s effect size metric (Lakens 2013). Hedges's g s is a bias corrected standardized median difference (Cohen's d s ) measuring the degree to which sample results depart from anticipations expressed in the null hypothesis being tested from the data (Thompson 2007). The following Hedges's g s thresholds for stating effect sizes were adopted: small (g s = 0.2); medium (g s = 0.5) and large (g s = 0.8) (Thompson 2007, Lakens 2013. For heatwaves, the low number of events and their accumulation exclusively in summer months (June-August), required a different approach. A summer was considered to experience a heatwave if the HWMID value exceeded zero on any day during that summer. Thus, a maximum of one heatwave event could be defined for each year. Subsequently, every summer was assigned with values for the water quality variables. Assuming that a heatwave would have a long-lasting imprint on water quality during the summer period, the value of water quality variable in August was assigned. Thus, each summer either had or did not have a heatwave, and was paired with epilimnetic and hypolimnetic water quality data in August, after which, the analyses proceeded in the same way as the analyses examining the impacts of extreme drought and wet events.
All calculations and graphical illustrations were implemented in (R Core Team 2017).

Historical changes in hydrometeorological and water quality variables
The period after WWTPs was characterized by lower medians of precipitation and inflow (Mann Whitney tests on pre-whitened series) than before WWTPs (figure S16, table S1). Interestingly, the frequency of extreme droughts, relative to total drought events, increased after WWTPs, from 5%-10% to 15%-16%, depending on the index (table S3(a)). Summers with heat waves also increased from 35% to 58% after WWTPs (table S3(c)). In contrast, the frequency of extreme wet events decreased after WWTPs, from 9%-10% to 4%-5%, depending on the index (table S3(b)). Before WWTPs, the median of hypolimnetic DO was 105% lower than its corresponding median after WWTPs. On the other hand, the median of NH 4 + was 98% higher than its corresponding median after WWTPs. Medians of NO 3 − , TP and T w were statisticslly similar between the two periods ( figure 1, table S2). In the epilimnion, the medians of NO 3 − and NH 4 + were 75% and 98%, respectively, lower while the median of TP was 35% higher after WWTPs, compared with their respective medians before WWTPs ( figure 1, table S2). Overall, the building and upgrading of WWTPs drastically reduced concentrations of nitrogen compounds in the reservoir, enhanced dissolved oxygen in the hypolimnetic water layers but did not seem to affect total phosphorus.

Causal relationships between hydrometeorological drivers and water quality variables
The building and upgrading of WWTPs affected the strength of Granger causality between hydrometeorology and water quality variables (figure 2, supplementary figures S3 and S4). The causal influence of inflow on hypolimnetic DO was particularly strong before WWTPs but disappeared after WWTPs (figures 2 and S3). In contrast, Granger causality of inflow on DO was altogether absent in the epilimnion in both periods (figures 2 and S4). After WWTPs, hypolimnetic NO 3 − showed causality with Tmax, covering the whole quantile range, and low to mid quantiles for epilimnetic NO 3 − and DO (figures 2, S3 and S4). TP was unresponsive to the four predictors, in both periods, except for a weak Granger causality with Tmax before WWTPs (figures 2 and S3). Unsurprisingly, Tmax showed a strong causality with T w , which became even stronger in the epilimnion after WWTPs (figures 2, S3 and S4).

Impact of ECEs on water quality before and after WWTPs
Overall, the hypolimnion manifested a clear reduction of the impacts of ECEs after WWTPs, whereas the response in the epilimnion defied simple description and suggested a complex set of responses both before and after WWTPs (figure 3, tables S4-S8).
In the hypolimnion, extreme (streamflow) droughts were associated with a lower median of DO (moderate Hedges's g s , figures 3 and S5, panel (C)) and higher medians for NH 4 + (moderate Hedges's g s , figures 3 and S7, panel (C)) and T w (weak Hedges's g s , figures 3 and S9, panel (C)) relative to non-extreme drought conditions, before WWTPs. However, after Figure 2. Normalized Granger causality of hydrometeorological predictors (air temp., inflow, precipitation, and water level) on water quality variables (DO = dissolved oxygen, NH4 + = ammonium, NO3 − = nitrate, TP = total phosphorus, and water temp = water temperature). The Granger causality analysis differentiates between the hypolimnion and epilimnion, and between the periods after and before WWTPs. For each analysis, the stacked bars show the normalized Granger causality integrated over three quantile regions (differentiated by three colours). Thus, for cases where no bars are present, there is no evidence of Granger causality for that combination of predictor, water quality variable, layer, and period. The presence of stacked bars imply that significant Granger causality was identified for the quantile region represented by the colour. For instance, air temperature showed Granger causality on NO3 − after WWTPs for both the epilimnion and the hypolimnion, and the causality accumulated in the quantile region between 0.25 and 0.75 (green bar). However, there is no evidence of causality between the same variables before WWTPs. For a detailed representation on how normalized Granger causality is calculated, see figures S2-S4.
WWTPs, none of the extreme drought conditions was associated with changes in any water quality variable. Extreme (inflow, water level & precipitation) wet events were associated with higher medians of DO (weak Hedges's g s , figures 3 and S5, panel (M)), NO 3 − (strong Hedges's g s , figures 3 and S6, panels (E), (M)) and T w (weak Hedges's g s , figures 3 and S9, panel (I)) and lower medians of NH 4 + (strong & moderate Hedges's g s, figures 3 and S7, panels (I), (M)) and TP (strong Hedges's g s, figure 3 and table S7), relative to their respective non-extreme wet conditions, before WWTPs. While similar (streamflow) events were associated with a higher median value of DO (strong Hedges's g s , figures 3 and S5, panel (F)), the rest of water quality variables, however, were unaffected by extreme wet conditions after WWTPs. Before WWTPs, summers with heatwaves were associated with a lower median of DO (weak Hedges's gs, figures 3 and S5, panel (A)) compared to summers without heatwaves, whereas the rest of the water quality variables were unresponsive to heatwaves. After WWTPs, summers with heatwaves were associated with a lower median of TP (strong Hedges's g s , figure 3 and table S7) compared to summers without heatwaves, whereas no other water quality variable responded to heatwaves (figure 3, tables S4-S8).
In the epilimnion, extreme (inflow & precipitation) drought events were associated with lower medians of DO (strong Hedges's g s , figures 3 and S10, panel (C)) and NO 3 − (strong Hedges's g s , figures 3 and S11, panel (K)) and higher medians of NH 4 + (weak Hedges's g s figures 3 and S12, panel (C)) and T w (strong Hedges's g s figures 3 and S14, panel (K)), relative to their respective non-extreme drought conditions, before WWTPs. Similar extreme (inflow & water level) droughts were not associated with any changes in DO but corresponded to a higher median of NH 4 + (weak Hedges's g s , figures 3 and S12, panel (D)) and a lower median of T w (strong Hedges's g s figure 3 and table S8) after WWTPs. Extreme (precipitation) wet events were associated with a higher median of NO 3 − (strong Hedges's g s , figures 3 and S11, panel (M)) relative to their nonextreme conditions before WWTPs. Similar extreme (precipitation) wet events were not associated with any change in NO 3 − but corresponded to a lower median of TP (strong Hedges's g s , figure 3 and table S7), after WWTPs. Furthermore, extreme (water level and inflow) wet conditions were associated with a lower median of NH 4 + (weak Hedges's g s figure 3 and table S6) and a higher median value of DO (strong Hedges's g s , figures 3 and S10, panel (F)), respectively, after WWTPs. Summers with heat waves were associated with a decrease in the median of NO 3 − (strong Hedges's g s , figures 3 and S11, panel (B)) relative to summers without heat waves after WWTPs, whereas no heat wave effect was observed on any water quality variable before WWTPs (tables S4-S8). After WWTPs, summers with heatwaves were associated with a lower median of TP (strong Hedges's g s , figure 3 and table S7) compared to summers without heatwaves, whereas no any other water quality variable responded to heatwaves (figure 3, tables S4-S8).

Discussion
The hypothesis that a change in trophic state of the reservoir increased resistance of the system's water quality to impacts of ECEs, was supported by results from both the Granger causality in quantile (particularly on DO) and behaviour of water quality under ECEs analyses. Both analyses also suggested that the change in resistance of water quality to ECEs, before and after WWTPs, was particularly prominent in the hypolimnion, whereas the epilimnion retained a comparable sensitivity to ECEs. The absence of statistically significant differences in medians of almost all water quality variables (except for DO), between extreme and non-extreme hydrological events, in the hypolimnion, i.e. the absence of impacts of hydrological extreme events on water quality, after WWTPs, is interpreted as increased resistance, following the definitions provided for by Perfecto et al (2019), Thayne et al (2021) and Mitra et al (2015. Conversely, in the epilimnion, the presence of statistically significant differences in medians of water quality variables between extreme and non-extreme hydrological events, in both periods, implies the retention of sensitivity to impacts of those ECEs, hence interpreted as lack of or less resistance. The increased resistance of water quality to ECEs in the hypolimnion, after WWTPs, coincided with an increased frequency of extreme droughts and heatwaves, making it unlikely that these observations were an artefact caused by a reduction in the number of ECEs. The differences in behaviour between the epilimnion and hypolimnion suggest that reservoir biogeochemical processes play a fundamental role in defining the response of water quality to ECEs. We postulate that the observed differences in Granger causality and impacts of ECEs on DO concentrations are governed by (a) differing metabolic pathways of anabolism and catabolism dominating the epilimnion and hypolimnion respectively, and (b) stronger water-air coupling in the epilimnion making it susceptible to transient hydrometeorological conditions. The dominance of primary production in the epilimnion may have influenced the observed lack of causality between hydrometeorological variables and water quality (except Tmax). Phytoplankton is a dynamic component of lake ecosystems, with a generation time in the order of days, whereas the monitoring program in Sau Reservoir was on a monthly basis. While a monthly sampling resolution can capture changes driven by eutrophication (Pomati et al 2012, Jochimsen et al 2013), it is insufficient to capture the short term dynamics of primary producers, operating on time scales of hours to days. Using monthly data may have limited the capability to predict water quality changes in the epilimnion.
Hypolimnetic water quality was more sensitive to ECEs probably due to the crucial role of organic matter degradation as a pathway for energy and matter flow in deep layers. Organic matter degradation in the hypolimnion is driven by the flow of organic materials from the river and the epilimnion, and the availability of electron acceptors (DO and NO 3 − ). Degradation of organic matter does not show short-term oscillations typical of primary production processes (Bastviken et al 2004), which may explain our observed causality between inflow and DO.
High causality and the significant impacts of ECEs found in the hypolimnion, before WWTPs, could also be explained by background water quality conditions. Before the operationalization of the WWTPs, the hypolimnion was anoxic, with a very high NH 4 + to NO 3 − molar ratio (Marcé et al 2008b). This scenario suggests a stronger role of organic matter degradation processes that keeps DO levels low. The high NH 4 + to NO 3 − ratio originated from the fact that most dissolved nitrogen entered the reservoir as NH 4 + and the presence of anoxia promoted denitrification (Hedin et al 1998, Burgin et al 2011. Thus, extreme drought conditions would exacerbate respiration processes, hence decreasing DO and increasing NH 4 + concentrations. In contrast, extreme wet conditions would favour re-oxygenation, increase the NO 3 − load from non-point sources, and dilute NH 4 + concentrations from WWTPs. After WWTPs, only re-oxygenation during extreme wet events remained, suggesting a reduced role of respiration processes in the hypolimnion, which was no longer responding to extreme droughts through low DO and high NH 4 + levels. Thus, the dramatic decrease in NH 4 + concentration blurred the dilution effect of extreme wet events.

Conclusion
This study provides a novel approach to infer causality between hydrometeorological extreme events and reservoir water quality using low-frequency water quality time-series. The non-parametric causalityin-quantile approach to analysing a long term, low frequency dataset, gathered from a real and complex system, shows that trophic state modulates the resistance of reservoir water quality variables to effects of ECEs and that the deeper water layers are more strongly impacted by these events and hence their resistance is more dependent on trophic state. We hypothesise that by keeping reservoirs in good trophic state (mesotrophic to oligotrophic), their water quality becomes more resistant to impacts of extreme events, giving water resources managers some control in dealing with adverse effects of the climate extremes. We conclude that besides improving water quality (in order to achieve UN Sustainable Development Goal 6.3), having water uncompromised by nutrient pollution also ameliorates adverse effects of climate extremes on reservoir water quality.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.