Long-term river invertebrate community responses to groundwater and surface water management operations.

River flow regimes have been transformed by groundwater and surface water management operations globally, prompting widespread ecological responses. Yet, empirical evidence quantifying the simultaneous effects of groundwater and surface water management operations on freshwater ecosystems remains limited. This study combines a multi-decadal freshwater invertebrate dataset (1995-2016) with groundwater model outputs simulating the effects of different anthropogenic flow alterations (e.g. groundwater abstraction, effluent water returns) and river discharges. A suite of flow alteration- and flow-ecology relationships were modelled that tested different invertebrate community responses (taxonomic, functional, flow response guilds, individual taxa). Most flow alteration-ecology relationships were not statistically significant, highlighting the absence of consistent, detectable ecological responses to long-term water management operations. A small number of significant statistical models provided insights into how flow alterations transformed specific ecological assets; including Ephemeroptera, Plecoptera and Trichoptera taxa which are rheophilic in nature being positively associated with groundwater abstraction effects reducing river discharges by 0-15%. This represents a key finding from a water resource management operation perspective given that such flow alteration conditions were observed on average in over two-thirds of the study sites examined. In a small number of instances, specific invertebrate responses displayed relative declines associated with the most severe groundwater abstraction effects and artificial hydrological inputs (predominantly effluent water returns). The strongest flow-ecology relationships were recorded during spring months, when invertebrate communities were most responsive to antecedent minimum and maximum discharges, and average flow conditions in the preceding summer months. Results from this study provide new evidence indicating how groundwater and surface water resources can be managed to conserve riverine ecological assets. Moreover, the ensemble of flow alteration- and flow-ecology relationships established in this study could be used to guide environmental flow strategies. Such findings are of global importance given that future climatic change and rising societal water demands are likely to further transform river flow regimes and threaten freshwater ecosystems.


Introduction
The biodiversity of freshwater ecosystems has declined rapidly in recent years, with such losses occurring at over double the rate of that experienced in terrestrial or marine environments ( Tickner et al., 2020 ). The management of freshwater resources in riverine environments through the construction of infrastructure, groundwater and surface water abstraction practices, and artificial flow releases is a primary driver transforming freshwater ecosystems globally ( Arthington et al., 2018 ). Hydrological alterations within fluvial environments are ubiquitous worldwide ( Wada et al., 2013 ), with just 37% of the world's longest rivers flowing without the presence of human infrastructure (Grill et al., 2019). As such, managing freshwater resources to balance societal and ecosystem needs is a key challenge of the 21 st century ( Naiman and Dudgeon, 2011 ;Arthington et al., 2018 ).
Biota inhabiting riverine environments have developed a suite of morphological, behavioural and life-history strategies over evolutionary timescales which allow them to inhabit rivers conveying naturally varying flow regimes ( Lytle and Poff, 2004 ). The transformation of river flow regimes threatens biota devoid of the necessary traits required to adapt to anthropogenic flow alterations ( Ruhi et al., 2018 ). As such, there have been numerous attempts to rehabilitate altered flow regimes and mitigate aquatic biodiversity losses globally ( Gillespie et al., 2015 ;Poff et al., 2017 ). However, riverine ecosystem responses to flow regime rehabilitation efforts have been unpredictable, inconsistent between studied taxa and highly variable across space and time ( Gillespie et al., 2015 ). Establishing ecological responses to hydrological controls (i.e. flowecology relationships) across large spatial and temporal scales is vital for guiding environmental flow (e-flow) methodologies aiming to balance societal and ecosystem water demands ( Davies et al., 2014 ). Various e-flow frameworks are guided by flow-ecology relationships that help identify specific hydrological events benefitting specific ecological assets (e.g. target taxa), such as the 'Functional Flows' ( Yarnell et al., 2020 ) and 'Designer Flow' methodologies ( Chen and Olden, 2017 ).
Studies quantifying ecological responses to anthropogenic flow alterations (i.e. flow alteration-ecology relationships) are more limited globally compared to those establishing flow-ecology relationships ( Bradley et al., 2017 ), which hinders our scientific understanding of how existing water resource management operations have influenced riverine ecosystems. Moreover, studies developing flow alteration-ecology relationships based on groundwater abstraction effects have been notoriously understudied ( Poff and Zimmerman, 2010 ). Instead, riverine ecological responses to groundwater abstraction are often inferred based on experimental designs manipulating water volumes associated with hypothetical low-flow/abstraction scenarios ( Dewson et al., 2007 ;Aspin et al., 2019 ). As a result, there remains a paucity of evidence quantifying ecosystem responses to empirically derived groundwater abstraction pressures ( Gleeson and Richter, 2018 ). Such ecological evidence and appraisals are urgently required given the pervasive effects of groundwater abstraction on riverine ecosystems, which have reduced river discharges by an average of 10% worldwide ( de Graaf et al., 2014 ), and greater societal water demands are projected to increase groundwater abstraction pressures in the future ( van Loon et al. 2016 ).
To our knowledge, this study provides the first evidence of long-term  ecological responses to anthropogenic flow alterations driven by both groundwater and surface water management operations, as well as historic hydrological (discharge) variations. To achieve this, freshwater invertebrate biomonitoring data was combined with regional groundwater model outputs to construct an ensemble of flow alteration-and flow-ecology relationships.

Study region
The study region is situated in the south-west of the United Kingdom (UK; Fig. 1 ) within the operational boundaries of the water company 'Wessex Water' plc. (WW), who manage water resources for 2.8 million inhabitants. The focus of this study was centred on groundwater dominated river systems predominantly underlain by Cretaceous Chalk (a white, fine-grained limestone -calcium carbonate). The underlying aquifer supports more groundwater abstraction than any other aquifer in the UK and accounts for 60% of the groundwater used across England and Wales ( Macdonald and Allen, 2011 ). In the WW region, 75% of water abstracted by is derived from groundwater sources across the study region (from 72 public water sources -Wessex Water, 2014 ). However, outflows from 97 effluent water returns results in some river reaches conveying greater flows than would occur naturally. In addition, the study region contains 6 low flow alleviation schemes, where groundwater is augmented into select rivers which fall below a designated discharge. The study region is dominated by arable agriculture and grassland ( National River Flow, 2020 ) and the rivers are typically characterized by high dissolved oxygen levels, alkalinity and nutrient levels ( White et al., 2018 ).

Wessex Basin groundwater model
The 'Wessex Basin' groundwater model (see Soley et al., 2012 ) was used to characterize anthropogenic flow alterations and hydrological variations at each study site. The Wessex Basin model divides the WW region underlain by Cretaceous Chalk and Upper Greensand into 250 × 250m grid cells. The Wessex Basin model has been adapted from the finite difference MODFLOW model (see McDonald and Harraugh, 1988 ) and estimates the interaction between stream cells and groundwater levels at circa ( c .) 10-day intervals (3 modelled outputs per month). These MODFLOW principles were then combined with daily outputs from a 4R (Rainfall, Recharge and Runoff Routing) hydrological model to estimate the total daily discharge ('historic' discharges herein) at each stream cell (see Heathcote et al., 2004 ). Historic discharge estimates incorporate groundwater abstraction influences (from WW and other small abstractors) and artificial hydrological inputs (principally effluent water returns, but also some low-flow alleviation strategies). The Wessex Basin model also models daily 'naturalised' discharges -those not subjected to any form of hydrological alteration. Anthropogenic flow alterations were derived from the percentage difference between the naturalised and historic discharge time series.
Anthropogenic flow alteration and historic discharge time series were obtained from 01/01/1994 to 31/12/2016 for each stream cell representing sampling sites (thus incorporating flow data from the 12-month period preceding all invertebrate samples; see below). Some short-term drying events identified at 9 study sites were incorrectly modelled according to measurements from nearby flow gauges; but zero-flow days totalled < 6% of the total hydrological time series at any given site and were subsequently replaced by interpolated values derived from the na.approx function in the 'zoo' package ( Zeileis, 2020 ). This was performed within R studio ( R Development Core Team, 2014 ) along with all subsequent statistical analyses.

Invertebrate dataset
Data on instream invertebrate communities was obtained from the 'BIOSYS' database collated by the Environment Agency (EAthe environment regulator for England). Invertebrate samples collected during WW's routine biomonitoring were also used in the analysis. Surveyors from both organizations adopted a standardized kick-sampling procedure, with all habitats within a river reach being proportionally sampled over a 3-minute period, supplemented with an additional 1-minute hand search ( ISO, 2012 ). Both EA and WW invertebrate samples were subjected to industrial standard The location of invertebrate sampling locations (circles) along rivers (black lines) in relation to hydrogeological conditions and their associated flow alteration (%) averaged across the study period. For hydrogeological conditions: dark grey = highly productive aquifer; light grey = moderately productive aquifer; white = low productivity aquifer or rocks with essentially no groundwater (for classification, see British Geological Survey, 2020 ). quality control assurances. Samples collected between 1995-2016 were utilized in the study.
Only study sites exhibiting a groundwater dominated hydrology, perennial flows (guided by the groundwater model and gauged flows) and accurate groundwater model outputs (guided by validation exercises of the groundwater model -ENTEC, 2015 ) were retained. In total, 1693 invertebrate samples collected during spring (March-May) and autumn (September-November) from 89 study sites containing at least 5 years of invertebrate data were initially selected, and region-wide hydrological information from these were analysed throughout. However, for ecological analyses, due to subsequent data refinement (removal of samples containing outliers -see Section 3.2 ), 1405 invertebrate samples (705 in spring and 700 in autumn) from 87 study sites were analysed, with a mean average of 19.2 samples and 10.1 years of data per study site.
Abundances of all invertebrate taxa were aggregated into logarithmic groups, which were subsequently recorded on an ordinal scale (so that 1 = 1-9 individuals; 2 = 10-99; 3 = 100-999; 4 = 10 0 0-9999; 5 = ≥10 0 0 0; sensu Durance and Ormerod, 2009 ). Invertebrates were predominantly identified to family-level in accordance with the taxonomic resolution routinely utilized by the EA ( ISO, 2012 ). Due to variations in how taxa were resolved over time, the following community data harmonization protocols were implemented: Cragoncytidae was combined with Gammaridae; Limonidae, Pedicidae and Scertidae were grouped under Tipulidae; Lumbricilidae, Naididae, Tubifidae, Haplotaxidae, Lumbricidae, Nematoda and Nematomorpha were pooled with Oligochaetes; Dugesidae were included in Planaridae, and Bithynidae were incorporated with Hydrobiidae. While coarser taxonomic resolutions may lead to weakened hydrological signatures on invertebrate communities ( Monk et al., 2012 ), Durance and Ormerod (2009) highlighted that invertebrate data resolved to the same taxonomic resolution yielded highly significant responses to antecedent discharges variations in the same region. Finally, Bibionidae, Collembola, Hebridae, Oribatei, Staphylinidae, Succineidae, Zonitidae were excluded from the analysis as they are typically indicative of terrestrial/riparian %EPT characterizes the percentage of taxa recorded in a sample belonging to Ephemeroptera, Plecoptera and Trichoptera orders.

Flow response guild
Family-level Loticinvertebrate Index for Flow Evaluation -'Family LIFE' The LIFE score is a biomonitoring index summarising the velocity preferences of invertebrate taxa ( Extence et al., 1999 ). LIFE is routinely used by UK regulatory bodies to set water abstraction licence conditions and monitor invertebrate responses to flow variability ( Klaar et al., 2014 ). Slower flow taxa richness comprise taxa in flow groups 3 and 4 (preferring slow flow velocities to standing waters) and Rheophilic taxa richness consist of taxa in flow groups 1 and 2 (preferring moderate to rapid flow velocities, > 20cm/s).
Richness of taxa preferring slow velocities -'Slower flow taxa richness' Richness of rheophilic taxa-'Rheophilic taxa richness' Functional Functional diversity Functional diversity was calculated using inverse Simpson's index performed on a trait-abundance array.

Rao's Quadratic Entropy
Rao's Quadratic Entropy (the product of the pairwise functional distances between taxa); Functional Richness (the minimum convex hull in functional space encompassing all species); Functional Evenness (the regularity of distances in functional space between taxa connected by a minimum spanning tree); Functional Divergence (the distribution of taxa in relation to the functional centroid -see Villéger et al., 2008 ) were calculated using the dbFD function in the 'FD' package ( Laliberté et al., 2015 ) and derived from a Bray-Curtis dissimilarity matrix.

Functional Richness Functional Evenness Functional Divergence
Individual taxa Ancylidae Although part of the Planorbidae family, this taxon is routinely recorded separately. Only 1 species comprises Ancylidae ( Ancylus fluviatilis ) in the UK and is rheophilic. Asellidae In the UK, 4 Asellidae species have been recorded, 2 of which are dominant in fluvial environments across the study region: Asellus aquaticus and Proasellus meriadanus . Both species prefer slow flow velocities or standing waters. Ephemerellidae Only 2 Ephemerellidae species have been recorded in the UK and both prefer moderate to fast flow velocities. In the study region, 1 species ( Serratella ignita ) widely occurs in fluvial environments. Ephemeridae In the UK, 3 Ephemeridae species have been recorded and all prefer moderate to fast flow velocities. Only 1 species is dominant in fluvial environments across the study region ( Ephemera danica ). Erpobdellidae In total, 5 Erpobdellidae species have been recorded in the UK and all prefer slow flow velocities or standing waters (except for Trocheta bykowskii , never sampled in the study region). Only 1 species ( Erpobdella octoculata ) widely occurs in fluvial environments across the study region. Glossiphoniidae In the UK, 8 Glossiphoniidae species have been recorded and all species prefer slow flow velocities or standing waters; with 3 species being commonly sampled in fluvial environments across the study region ( Theromyzon tessulatum, Glossiphonia complanata and Helobdella stagnalis ).

Rhyacophilidae
In the UK, 4 Rhyacophilidae species have been recorded which all prefer rapid flow velocities.
Only 1 species widely occurs in fluvial environments across the study region ( Rhyacophila dorsalis ). Sericostomatidae Only 2 Sericostomatidae species have been recorded in the UK, with 1 being most common in fluvial environments across the study region ( Sericostoma personatum ) and prefers moderate to fast flow velocities.
N.B. The velocity preferences of individual taxa were derived from flow groups underpinning the Lotic-invertebrate Index for Flow Evaluation (LIFE) score (see Extence et al., 1999 ). The occurrence of species in the study region are based on records in the National Biodiversity Network (NBN -https://nbn.org.uk/ ), existing regional databases (limited species-level samples in BIOSYS and Wessex Water invertebrate databases) and other reports and published literature ( Pardo and Armitage, 1997 ;Wessex Water, 2008 ;Armitage and Bass, 2013 ;White et al., 2018 ;White et al., 2019b ; Wessex Chalk Streams and Rivers Trust. 2020 ).
habitats. Following harmonization, 102 taxa were used in the analyses.

Invertebrate community responses
In total, 19 ecological response metrics were examined that characterized various community responses (characterizing taxonomic, functional and flow response guild properties) and individual taxa abundances (see Table 1 ).
The functional diversity of invertebrate communities was derived from the European database compiled by Tachet et al (2010) which adopts a fuzzy-coding procedure. Within this, faunal affinities to individual traits range from zero to three or five (indicating no to high affinity -the upper limit depending on existing scientific certainty). In total, 11 grouping features (i.e. a functional trait category -e.g. 'maximum body size') comprising 63 traits (i.e. modalities residing within grouping features -e.g. ' ≤0.25cm', ' ≥8cm') were examined that characterize the biological properties of invertebrate taxa (see Supplementary Material, Appendix A). Taxa that did not occur within the UK (guided by Davies and Edwards, 2011 ) were removed from the database, along with Chironomidae and all taxa resolved beyond family-level. Trait values were then standardized for all samples so that each grouping feature summed to 1 (to ensure trait affinities had equal weighting between taxa). Family-averaged trait values were then calculated and standardized (as above) to account for taxa expressing no affinity for all traits within a specific grouping feature. These standardized values were used to derive all functional metrics except for functional trait diversity (see Table 1 ), which was calculated from a trait × abundance array. For this, standardized, family-averaged trait values were multiplied by community abundances, averaged across all sampled taxa and standardized once more across all grouping features.
The individual taxa examined comprise UK species representatives expressing the same ecological guilds (defined by their preferences towards flow velocities -sensu Extence et al., 1999 ) and occurred widely throughout the invertebrate dataset (68%-79% of samples). In summary, 19 ecological responses were tested, 11 community response metrics and 8 individual taxa ( Table 1 ).

Hydrological and anthropogenic flow alteration temporal variations
Temporal variations in anthropogenic flow alterations and historic discharges were explored using values averaged on a monthly basis across all 89 groundwater stream cells (corresponding to invertebrate sampling locations). For anthropogenic flow alteration time series, a Generalized Additive Model was constructed using the sample date as a smoothing term so the degree of smoothness equalled 0.3 times the number of years using the gam function in the 'mgcv' package ( Wood, 2018 ) and all residual plots were inspected to ensure the assumptions of normality and homoscedasticity were satisfied.
For historic discharges, a flow regime magnitude (RM) classification procedure was undertaken to characterize periods of above and below average discharge conditions. Daily historic discharge values from each stream cell were transformed to z-scores. For the RM classification, historic discharge time-series were divided into hydrological years commencing in August and terminating in July to help ensure that the rising limb, annual peak and flow recession were typically incorporated across a 12-month period ( sensu Monk et al., 2006 ). Subsequently, monthly summaries of historic discharges (mean, standard deviation, minimum and maximum values) were obtained for each hydrological year. These 4 statistical summaries underwent unit-based standardization (X' = (X-X min /X max -X min ) + 1) to ensure equal weighting and were inputted into a hierarchical, agglomerative cluster analysis (Ward's method; see Monk et al., 2006 ). Subsequently, hydrological years were grouped into 4 flow magnitude classes (guided by Silhouette analysis implemented via the fviz_silhouette function in the 'factoextra' package - Kassambra and Mundt, 2017 ) spanning from low (RM1), moderate (RM2 and RM3) and high flow conditions (RM4).
Finally, to test the association between groundwater and surface water management operations and river flow regimes, anthropogenic flow alteration and historic discharge (z-score) values were averaged on a monthly basis and across all 89 stream cells; which were subsequently used as dependent and independent variables (respectively) in a linear regression using a logarithmic function. For this, residual diagnostics were inspected to ensure the data was normally distributed and homoscedasticity existed.

Invertebrate responses to hydrological and flow alteration indices
A suite of hydrological indices characterizing different river flow regime components were calculated, including 33 hydrological indices outlined in the 'Indicators of Hydrological Alteration' methodology ( sensu Richter et al., 1996 ) and additional ecologically important indices identified in previous research in the UK ( Worrall et al., 2014 ;White et al., 2019b ). These indices were calculated for both anthropogenic flow alteration ('AF' -n = 43; only AF indices characterizing flow magnitude alterations were processed) and historic discharge ('Q' -n = 47) time series (see Supplementary Material, Appendix B). Hydrological indices with skewed distributions or a lack of unique values ( < 100) were excluded from the analyses (n = 12). Subsequently, all hydrological indices underwent a unit-based standardization (X' = (X-X min /X max -X min ) + 1) and extreme outliers were identified from interquartile range (IQR) values (values that fall below Q1 − 3 ×IQR or above Q3 + 3 ×IQR). Subsequently, 10 indices were removed that contained > 50 outliers from the analysis (see Supplementary Material, Appendix B) and any sample containing an outlier within the remaining AF and Q indices were not analysed, which resulted in 1405 samples (and 87 sites with ≥5 years of data) being analysed.
In total, four statistical model groups were constructed that explored the separate influence of AF (i.e. flow alteration-ecology relationships) and Q (i.e. flow-ecology relationships) indices on invertebrate community responses derived from spring and autumn samples (AF-Spring; AF-Autumn; Q-Spring and Q-Autumn; see Fig. 2 ). Each of these model groups contained a series of Quantile Regressions (QRs) and Quantile Mixed-Effect Regressions (QMRs), which were performed across various quantiles (from 0.05 to 0.95 in 0.05 increments) to test community responses across a range of data values (including central tendencies, and upper or lower limits). QMRs incorporated river identity as a random effect to account for potential spatial autocorrelation and that samples from the same watercourse likely being correlated over time ( Mathers et al., 2020 ). QRs and QRMs were implemented using the lqm and lqmm functions (respectively) within the 'lqmm' package ( Geraci, 2020 ).
For each community response metric, a set of QRs and QRMs were initially performed for each individual hydrological index (independent variable; see Fig. 2 ). Each of these sets comprised the same 6 null models -QRs and QRMs that tested a constant response, spatial and/or temporal effects (i.e. no hydrological index) -alongside 26 'Single Hydrological Models' (SHMs) that modelled the influence of a single hydrological index via 4 statistical functions (Linear, Exponential, Logarithmic and Quadratic -sensu Fornaroli et al., 2019 ) together with different spatial and temporal effects (all formulae underpinning QRs and QRMs performed for each set are displayed in Supplementary Material, Appendix C). Akaike weights ( w i -derived from corrected Akaike Information Criteria values) were calculated and averaged across all quantiles for all QRs and QRMs ( Allen and Vaughn, 2010 ). A hydrological index would be excluded from subsequent analyses if its optimal (i.e. highest average w i ) SHM exhibited an average w i < 2 times the average w i of the optimal null model. The construction of SHMs allowed ecologically important hydrological indices, in the form of their most powerful statistical function, to be identified prior to multiple regression analyses (see below). Qualifying hydrological indices for each community response metric were then checked for collinearity by iteratively removing hydrological indices containing the highest variance inflation factor values until all were below 3.
Various 'full' Multiple Hydrological Models (MHMs) were created for each community response metric that modelled the additive effects of qualifying hydrological indices (modelled via their optimal statistical function identified in SHMs) alongside different combinations of spatial and temporal effects (the same as inputted into SHMs; see Fig. 2 ). Of these, the full MHM yielding the best variance structure (i.e. the model containing the spatial and temporal effects exhibiting the highest average w i ) was identified; which was subsequently compared against restricted models generated by removing AF or Q indices in a backward-stepwise procedure. This process was repeated iteratively until the 'optimal' MHM exhibiting the highest average w i was identified. Finally, the optimal MHM obtained for each community response metric was compared against a null model exhibiting the same spatial and temporal effects, but devoid of any Q or AF indices in order to test its overall significance (a method analogous to a likelihood ratio test). For this, significant, optimal MHMs (soMHMs) were identified as those yielding a relative average w i > 0.75 ( sensu Fornaroli et al., 2019 ).
For counts of individual taxa, the statistical formulae and procedures utilised in the QR and QRM approach described above were replicated using quantile count models, which were developed in accordance with Cade and Dong (2008) . The counts of invertebrate taxa were transformed by adding a random uniform number in the interval 0-1, (Y' = (Y + U[0,1)). Each QR and QRM was modelled across 100 different datasets comprising randomly added numbers, with the estimated coefficients being averaged to remove the source of additional variation introduced by adding random U[0, 1) numbers to Y.

Long-term hydrological variability and extreme flow conditions
Rivers in the study region displayed consistent intra-annual changes in hydrological variability (i.e. historic discharges from the groundwater model), with peak flows occurring during late winter to early spring, before declining across the summer and autumn months ( Fig. 3 a). Distinct periods of low (1996-1997, 20 04-20 06, 2010-2012) and high (20 01-20 01, 2012-2014) discharges were observed across the study region, which corresponded with the lowest (RM1) and highest (RM4) flow regime magnitude classes ( Figs. 3 a and 3 b). When anthropogenic flow alterations were averaged across all 89 study sites on a monthly basis, values ranged between -15-0% ( Fig. 3 a), and 67% of the study sites averaged flow alterations within this range across the entire study period. Anthropogenic flow alterations displayed the most extreme negative values at the beginning of the study period, but from 1997 onwards displayed inter-annual variations congruent with historic discharges, with the greatest reductions in discharge occurring during periods of low-flow (and vice versa ; Fig. 3 a). A highly significant logarithmic relationship was observed between the percentage of flow altered versus historic discharges (linear regressionr 2 = 0.59, F = 395, p < 0.001; see Fig. 3 c).

Flow alteration-and flow-ecology results structure
The remaining results highlight flow alteration-and flowecology relationships derived from the four statistical model 'groups' containing data from different seasons (spring and autumn) and hydrological index types (i.e. anthropogenic flow alteration -AF and discharge variability -Q). Each of these groups tested the effects of hydrological indices on each individual ecological response metric, therefore each group potentially contained up to 19 significant, optimal Multiple Hydrological Models (soMHMs;    the final step in the analytical framework adopted in this studysee Fig. 2 ).

Flow alteration-ecology relationships
Of the four model groups tested, AF-Spring and AF-Autumn (i.e. the two groups comprising flow alteration-ecology relationships) yielded the greatest number of non-significant models (each with 15 out of a possible 19; Table 2 ). Ecological response metrics yielding non-significant associations varied slightly between seasons, although the taxonomic diversity, 'Slower flow taxa richness' (the richness of taxa preferring slow-flow conditions), functional richness and functional evenness (see Table 1 for ecological response metric descriptions) and all individual taxa modelled (except for Ephemerellidae during autumn) were not significantly affected by AF indices ( Table 2 ).
The 4 significant, optimal Multiple Hydrological Models (soMHMs) in both the AF-Spring and AF-Autumn model groups comprised 8 and 9 AF indices, respectively. Of these, 10 AF indices characterized the discharge percentage altered during specific months, 7 of which yielded quadratic responses. In total, 5 of these modelled 'Rheophilic taxa richness' (the richness of rheophilic taxa) and '%EPT' (the percentage of Ephemeroptera, Plecoptera and Trichoptera taxa), which peaked at c. -15-0% flow alterations that characterise negligible-moderate abstraction intensities ( Figs. 4 a-d -58% to 87% of data values comprising these soMHMs lay between -15-0% AF values; see Supplementary Material, Appendix D, for figures displaying all associations comprising soMHMs not presented in the main text). However, it should be noted that the functional diversity and Rao Quadratic Entropy values displayed the opposite trends in response to September flow alterations based on invertebrate samples collected during spring and autumn, respectively (see Appendix D,Figures D1b and D2c). In total, 6 AF indices comprising soMHMs characterised the magnitude and duration of specific flow alterations, all of which yielded a logarithmic relationship. A key example of such a trend was 'Family LIFE' (the Family-level Lotic-invertebrate Index for Flow Evaluation) scores, which decreased most dramatically when the minimum flow alteration (i.e. most severe groundwater abstraction effect) in the 180-days prior to sampling (spring samples - Fig. 4 e) and 30-day average minimum flow alteration values (autumn samples - Fig. 4 f) dropped below c. -60% and -50% (respectively), which occurred in 1% of data values in both models.

Flow-ecology relationships
The Q-Spring model group yielded the highest number of soMHMs (n = 11), with the 9 non-significant models primarily comprising taxonomic responses and individual taxa ( Table 3 ). These soMHMs consisted of 32 Q indices, 14 of which characterised the duration of certain flow magnitudes, notably antecedent minimum and maximum discharges ( Table 3 ). One key soMHM association was Functional Evennessresponding unimodally to the maximum 7-day averaged discharge (Q7Max; Fig. 5 a). Some individual taxa responses also displayed important trends, such as Ephemerellidae abundances displaying a negative, linear trend with the maximum discharge in the 7-days prior to sampling (QMax7; Fig. 5 b) and Rhyacophilidae abundances increasing at faster rates with rising maximum 90-day average discharges (Q90Max; Fig. 5 c). In the Q-Spring group, 13 indices comprising soMHMs characterized the average discharges during specific months. Of these, 6 reflected summer flow conditions ( Table 3 ) and were influential on flow response guilds, with Family LIFE re-sponding in a positive, linear fashion to the average discharge in July (QJul; Fig. 5 d). Rheophilic taxa richness yielded a positive, logarithmic response to QJul ( Fig. 5 e), while Slower flow taxa richness displayed a quadratic response and notably declined across intermediate to high QJul values ( Fig. 5 f). The Q-Autumn group yielded 6 soMHMs that comprised 13 Q indices, which are summarised in Supplementary Material, Appendix D. Table 3 Optimal Multiple Hydrological Model outputs indicating the effects of historic discharge (Q) indices on ecological response metrics (i.e. flow-ecology relationships) derived from spring invertebrate communities. Flow regime components: M = Magnitude, F = Frequency, D = Duration, T = Timing and R = Rate of change. '1|' denotes a random effect in a quantile mixed-effect model. LIN = Linear, LOG = Logarithmic, QUA = Quadratic. Average wi denotes the mean Akaike weights across all studied quantiles. Significant, optimal Multiple Hydrological Models (soMHMs) are highlighted in bold.

Invertebrate responses to groundwater and surface water management operations
This study utilized a novel combination of long-term groundwater model outputs and a suite of invertebrate responses derived from regional biomonitoring programmes. While a limited number of studies have combined riverine ecological data with hydrolog-ical models characterizing groundwater abstraction effects across seasonal ( Kennen et al., 2014 ;White et al., 2019b ) to inter-annual periods ( Bradley et al., 2017 ;White et al., 2018 ;Liu et al., 2020 ), to our knowledge, the present study represents the first example to combine such sources of information over a multi-decadal temporal scale. Moreover, the groundwater model used accounted for different forms of water management operations, principally groundwater abstraction and effluent water returns (as well as a small number of low-flow alleviation strategies), which have seldom been explored simultaneously.
The present study highlighted that many flow alterationecology models were not statistically significant, including different taxonomic and functional community measures, as well as most individual taxa. Such findings highlight the absence of consistent, detectable ecological effects of long-term water management operations over a > 20-year period. This suggests that ecological responses to flow alterations were either highly variable (and therefore statistically undetectable) or exerted a weak influence on instream communities. On a national scale, Mims and Olden (2013) attributed weak statistical associations between daminduced flow alterations and fish communities to nuanced ecological responses that varied between different dam operations and flow regimes. However, the present study examined river systems at the regional-scale which yield comparable, groundwaterdominated flow regimes and were mostly subjected to lowmoderate groundwater abstraction intensities (approximately twothirds of study site averages and altered flow index values lay between -15-0%); as such, it is more likely that non-significant flow alteration-ecology relationships in the present study were an artefact of water resource management operations largely ex-erting weak (and not variable) ecological effects. Despite this, a small number of significant flow alteration-ecology relationships (8 in total -4 out of 19 in each season) provided some insight into how long-term groundwater and surface water management operations have modified certain ecological properties. For instance, the number of rheophilic taxa sampled during spring peaked when the preceding summer and winter discharges had been artificially reduced by c . 0-15%. This could be explained by groundwater abstraction practices in the winter immediately preceding the sampling period mitigating the effects of extreme highflow events on riverine biota (even rheophilic taxa; Miller et al., 2007 ), and those in the previous summer potentially concentrating flows and creating localized pockets of enhanced flow velocities ( White et al., 2019a ;2019b) . In addition, the percentage of Ephemeroptera, Plecoptera and Trichoptera (EPT) taxa sampled in autumn also peaked when groundwater abstraction effects reduced antecedent summer and winter discharges by c . 0-15%. Potential ecological mechanisms driving this may mirror those described for rheophilic taxa (as many representatives belong to EPT orders), in addition to low-moderate groundwater abstraction influences during the summer immediately preceding autumn samples potentially prompting earlier emergence and adult oviposition patterns, thus allowing new larval cohorts to develop in the autumn ( Brown et al., 2012 ). Equally, for both rheophilic and EPT taxa, affiliations with low to moderate groundwater abstraction effects could be attributed to greater invertebrate densities through organisms being concentrated into a reduced wetted perimeter ( Dewson et al., 2007 ), or channel contraction enhancing the sampling detection of a greater number of taxa. Whichever the underlying ecological mechanism(s) are, these unimodal flow alterationecology relationships provide vital empirical evidence highlighting that over a > 20-year time period, regional groundwater and surface water management operations broadly have not limited important riverine ecosystem components. However, some functional diversity measures displayed the opposite unimodal trend in relation to flow alterations during September, with relative declines in values at c . 0-15% groundwater abstraction effects potentially being due to greater dominance of EPT communities, which at the community-level yield more comparable functional niches ( Usseglio-Polatera et al., 20 0 0 ). Such findings have important implications globally, where there remains a paucity of information quantifying ecological responses to groundwater abstraction. For this reason, Gleeson and Richter (2018) proposed a global 'presumptive standard' on how groundwater resources should be managed to conserve surface water ecosystems and highlighted that baseflows (the portion of flow derived from groundwater or other delayed hydrological sources) should be reduced by no more than 10%.
While groundwater and surface water management operations in most instances reduced river discharges by c . 0-15%, some exceptions to this did occur and prompted distinct ecological responses. For instance, artificial hydrological inputs, the majority of which stemmed from effluent water returns, surprisingly led to relative declines in the number of rheophilic and EPT taxa sampled in spring and autumn, respectively (although the latter was less sensitive to flow inputs during winter). It is possible that alternative stressors associated with such water releases (e.g. modified water physio-chemical characteristics) were not suitable for sensitive biota (although it should be noted that major effluent water returns have been subjected to phosphorous stripping across the study region - Mainstone et al., 2005). In addition, some ecological response metrics, most notably EPT taxa (sampled during autumn) and Family LIFE, displayed relative declines when groundwater abstraction effects resulted in more severe negative flow alterations. While the occurrence of such extreme groundwater abstraction effects were rare ( c. 1% in the Family LIFE models), hydrological analyses in this study indicated that groundwater abstraction practices had a stronger relative effect on low-flow conditions; which occurs due to groundwater abstraction having a greater proportional influence on lower water levels that coincides with peaks in societal water demands. Given that combined increases in global temperatures and societal water demands are likely to exacerbate future low-flow conditions worldwide ( van Loon et al., 2016 ), further research is required to elucidate the direct effects of groundwater abstraction effects on riverine ecosystems.

Long-term structural and functional invertebrate responses to hydrological variations
The suite of flow-ecology relationships derived from such a dense network of ecohydrological information in this study has important implications for guiding water resource management operations and developing environmental flow strategies ( Davies et al., 2014 ). This study highlighted the importance of antecedent flow (discharge) variability for invertebrate communities sampled during the spring months. While previous studies within the UK have highlighted stronger flow-ecology relationships when examining invertebrate communities sampled during autumn (e.g. Monk et al., 2006 ;Worrall et al., 2014 ), Durance and Ormerod (2009) reported that antecedent river discharges yielded a greater influence on spring invertebrate communities in the same region as the present study. Such findings may reflect spring being a critical time period for the life-cycle of many invertebrates that are closely linked with antecedent flow conditions, including the maturation of insect larvae / nymphs prior to late spring and early summer emergence (e.g. Ephemera Danica -Everall et al., 2015 ) and periods of reproduction for others (Ancylus fluviatilis - Pfenninger et al., 2003 ).
For invertebrate communities sampled during spring, the duration of certain flow magnitudes (notably minimum and maximum discharges) were of critical importance. This includes functional evenness responding unimodally to the maximum 7-day averaged discharges, which potentially reflects competitively superior taxa dominating when antecedent high-flow conditions are less severe, while only highly tolerant taxa persisting when the maximum 7day averaged discharge is at its greatest (see Townsend et al., 1997 ). This has important implications for wider ecosystem functioning ( Villeger et al., 2008 ) and could inform environmental flow strategies that may target intermediate flow magnitudes to promote greater functional evenness of biotic communities. The abundance of Ephemerellidae during spring months declining with the maximum discharge in the 7-days prior to sampling could be due to their affiliation with fine-leaved macrophyte habitats, which become dominated by select taxa when hydraulic pressures intensify between plant strands (see White et al., 2019b ). This has important management implications given the socioeconomic importance of this mayfly for commercial and angling activities ( Everall et al., 2018 ). Conversely, the trichopteran Rhyacophilidae was positively associated with the maximum 90-day averaged discharges, emphasising the torrenticole nature of this key predator which benefits of long periods of sustained higher discharges ( Extence et al., 1999 ).
This study highlighted the sensitivity of specific invertebrate community flow response guilds to antecedent flow conditions, which supports the findings of Chen and Olden (2018), who highlighted the effectiveness of such responses in facilitating spatially transferrable flow-ecology relationships. From a national perspective, the sensitivity of Family LIFE during spring months represents an important finding as this biomonitoring index is used to guide water resource management operations and set abstraction licenses ( Klaar et al., 2014 ). Importantly, summer discharges had a strong effect on flow response guild values of spring invertebrate communities in the following year and highlights the long-lasting ecological effects of low-flow conditions ( Extence et al., 1999 ).

Conclusions
This study provides a novel perspective on the ecological effects of water resource management operations by combining biomonitoring datasets with groundwater model outputs over a multi-decadal timeframe. Our findings provide empirical evidence regarding how regional water resource management operations largely resulted in statistically non-significant effects on benthic invertebrate communities. However, a small number of statistical models recorded invertebrate communities being sensitive to more extreme groundwater abstraction intensities and artificial hydrological inputs, although such flow alterations were less common over the study period. Given a global paucity of evidence highlighting the long-term implications of groundwater abstraction effects on riverine ecosystems, studies like this are urgently required to inform the management of groundwater and surface water resource operations. Such research will become increasingly important in the face of increasing hydrological pressures from rising societal water demands and projected climatic changes.