Current and Future Ecological Status Assessment: A New Holistic Approach for Watershed Management

The Paiva River catchment, located in Portugal, integrates the Natura 2000 network of European Union nature protection areas. Resorting to topography, climate and land-use data, a semi-distributed hydrological model (Hydrological Simulation Program–FORTRAN) was run in order to simulate the hydrological cycle of the river and its tributaries. The model was calibrated over a 25-year period and validated within a 31-year period. Its performance was verified by comparing the recorded and simulated daily flows. The values of the Nash–Sutcliffe coefficient of efficiency of 0.95 and 0.76, and coefficient of determination of 0.95 and 0.82, were achieved for calibration and validation, respectively, thus showing a quite satisfactory model performance. Subsequently, the climate change impacts on temperature and precipitation, as well as their extremes, and on the flowrates were also assessed for a future period (2041–2070) under two anthropogenic forcing scenarios (representative concentration pathways 4.5 and 8.5). A procedure for selecting the most relevant metrics for assessing the ecological condition of the Paiva River was developed based upon a set of 52 invertebrate families sampled. Correspondence analyses were carried out for biological datasets (traits/metrics) with physicochemical and land use/land cover matrices separately. Out of all variables, water quality and flow and agriculture land use explained most of the variance observed. The integrated analysis undertaken in the present study is an important advance when compared to previous studies and it provides key information to stakeholders and decision-makers, particularly when planning suitable adaptation measures to cope with changing climates in the forthcoming decades.


Introduction
Natural integrity of river basins and the potential impact of anthropogenic pressures on water resources is raising social awareness worldwide, particularly regarding their socioeconomic prominence [1][2][3]. Hydrological modelling and simulation aim to quantitatively answer questions regarding the behaviour of water within a given catchment basin and assess its ecological status. and human pollution". This idea was emphasised by Bonada et al. [51], which found that studies using biological functional traits (at larger spatial scales and higher biological organizational levels) may anticipate the impacts of climate change. Also, Charvet et al. [52], Phillips [53], and Vieira et al. [54] stressed that, unlike taxonomic composition, biological traits are relatively stable across large environmental gradients (e.g., geology, altitude, geographical coordinates, stream order, and slope), allowing for their use for biomonitoring over large geographical areas. Accordingly, the evaluation protocols must be based on both conventional attributes and trait-based variables once together are sufficiently sensitive to distinguish impairment. Undeniably, indices that are only based on traits or together with traditional taxonomic methods are the future of the assessment of causal relationships between specific stressors and macroinvertebrate community response [55,56].
This study assesses the current water and ecological status of the Paiva River catchment, in Portugal, and the corresponding potential modifications under climate change projections. In order to achieve this goal, the objectives of this study are three-fold: (1) to calibrate and validate flowrate and water quality in the Paiva River catchment for the historical period 1950-2005; (2) collect climatic datasets from a five-member ensemble of GCM-RCM chain simulations in order to assess the impact of climate projections in flowrate and water quality under two different anthropogenic forcing scenarios (Representative Concentration Pathways, RCP), RCP4.5, and RCP8.5; and, (3) to demonstrate the importance of the macroinvertebrate traits in discriminating the different levels of water ecological status and to identify which of them can be more effectively and accurately used in future monitoring of the Paiva River catchment under a climate change context. Lastly, some preventive measures that are devoted to the conservation of this sensible aquatic ecosystem will be presented contributing future management measures.

Study Area
The Paiva River (Portugal) is part of Natura 2000 network (code PTCON0059), of nature protection areas in the territory of the European Union. It flows between two other Sites of Community Importance: "Serra de Montemuro" (to the north) and "Serras da Freita e Arada" (to the south) ( Figure 1). It is one of the Douro River main tributaries of the left margin and covers an area of approximately 796 km 2 , with a total extension of nearly 115 km [57] (Figure 1). The Paiva River catchment shows a very complex orography, with the spring at an elevation of approximately 1000 m and the river mouth at 10 m altitude, flowing in a deep canyon, and nearly halfway between very steep mountains. The land-use is mainly comprised of shrubs (38%), plantations of Pinus Pinaster (17%) and Eucalyptus (16%), followed by agriculture (14%). In its final section, the slopes show high coverage and good plant density despite the increase in the area planted with eucalyptus, already reflecting humid maritime conditions [58]. Temperate Mediterranean climate characterizes the catchment, with an average annual temperature of 13 • C and average annual precipitation greater than 1000 mm [45,46,59].

Climate Data
Very high resolution (~1 km) gridded daily mean, minimum, and maximum temperatures over mainland Portugal (PT.TG.HRES, PT.TN.HRES and PT.TX.HRES) over the baseline period of 1950-2015 were retrieved for the grid points covering the Paiva River catchment. The methodology for the development of these datasets is based on the E-OBS observational gridded dataset [24], as the underlying station data were already subjected to quality control procedures and homogeneity analysis. A two-step approach was delineated: (1) geostatistical interpolation of temperatures to capture daily mean spatial variability and (2) bilinear interpolation of daily temperature anomalies to capture temporal variability. The complete methodology followed to obtain these gridded datasets, their validation and characterization can be found in Fonseca and Santos [59]. Gridded daily precipitation totals at~1 km spatial resolution (PT.P.HRES) over the same baseline period were also retrieved for the Paiva River catchment. The PT.PHRES was developed based on gridded daily precipitation totals at~20 km spatial resolution (PT02 database) [60], supplied by the Portuguese Weather Service (www.ipma.pt). The steps that were followed for the calculation of the gridded dataset PT.P.HRES can be found in Fonseca and Santos [21]. These very high-resolution datasets provide detailed information on temperature and precipitation spatial and temporal variability in the Paiva River catchment, which will be an important added value for the subsequent hydrological simulations.

Climate Data
Very high resolution (~1 km) gridded daily mean, minimum, and maximum temperatures over mainland Portugal (PT.TG.HRES, PT.TN.HRES and PT.TX.HRES) over the baseline period of 1950-2015 were retrieved for the grid points covering the Paiva River catchment. The methodology for the development of these datasets is based on the E-OBS observational gridded dataset [24], as the underlying station data were already subjected to quality control procedures and homogeneity analysis. A two-step approach was delineated: (1) geostatistical interpolation of temperatures to capture daily mean spatial variability and (2) bilinear interpolation of daily temperature anomalies to capture temporal variability. The complete methodology followed to obtain these gridded datasets, their validation and characterization can be found in Fonseca and Santos [59]. Gridded daily precipitation totals at ~1 km spatial resolution (PT.P.HRES) over the same baseline period were also retrieved for the Paiva River catchment. The PT.PHRES was developed based on gridded daily precipitation totals at ~20 km spatial resolution (PT02 database) [60], supplied by the Portuguese Weather Service (www.ipma.pt). The steps that were followed for the calculation of the gridded dataset PT.P.HRES can be found in Fonseca and Santos [21]. These very high-resolution datasets provide detailed information on temperature and precipitation spatial and temporal variability in the Paiva River catchment, which will be an important added value for the subsequent hydrological simulations.
Gridded climatic datasets from a five-member ensemble of GCM-RCM chain simulations (see Table S1, Supplementary Material), developed within the framework of the EURO-CORDEX project, were retrieved for the development of climate change projections (http://www.euro-cordex.net/). The data were also subject to a bias correction performed under the previous project (SMHI-DBS45-MESAN, 1989-2010). The original spatial resolution is of about 12 km (EUR-11, 0.11° latitude × longitude) on a Gaussian mesh. The gridded daily precipitation totals (RR), minimum (TN), and maximum (TX) 2-meter air temperatures were extracted for the future period of 2041-2070 over a geographical sector covering the Paiva River catchment and under two different scenarios, RCP4.5 and RCP8.5. The data were then bilinearly interpolated to a latitude-longitude regular grid. Gridded climatic datasets from a five-member ensemble of GCM-RCM chain simulations (see Table S1, Supplementary Material), developed within the framework of the EURO-CORDEX project, were retrieved for the development of climate change projections (http://www.euro-cordex.net/). The data were also subject to a bias correction performed under the previous project (SMHI-DBS45- MESAN, 1989MESAN, -2010. The original spatial resolution is of about 12 km (EUR-11, 0.11 • latitude × longitude) on a Gaussian mesh. The gridded daily precipitation totals (RR), minimum (TN), and maximum (TX) 2-meter air temperatures were extracted for the future period of 2041-2070 over a geographical sector covering the Paiva River catchment and under two different scenarios, RCP4.5 and RCP8.5. The data were then bilinearly interpolated to a latitude-longitude regular grid.
Several steps were undertaken for downscaling the climate model data to~1 km spatial resolution: (1) for each model, scenario, and variable, data were bilinearly interpolated from~12 km to~1 km grid; (2) the gridded differences of monthly means (monthly sum) of temperature (precipitation) between model and observed data (PT.TG.HRES, PT.TN.HRES, and PT.TX.HRES) were calculated for a common baseline (1981-2010); (3) the corresponding long-term monthly differences at each grid point were then linearly interpolated to the daily timescale and smoothed out by a low-pass Lanczos filter, with a cut-off frequency at 30 days and 500 coefficients, to remove high frequency noise; and, (4) the gridded long-term daily differences were subsequently added (multiplied) to the simulated temperature (precipitation) for the period 2041-2070, under RCP4.5 and RCP8.5.

Hydrological Simulations
The climatic data were subsequently used as input in the Hydrological Simulation Program FORTRAN (HSPF) model in order to assess daily flowrates in the Paiva River. A graphical user interface of HSPF is included in the Better Assessment Science Integrating Point and Non-point Sources (BASINS), which is designed to perform water quantity and quality-based studies and was developed by the United States Environmental Protection Agency (USEPA) (https://www.epa.gov/ceam/betterassessment-science-integrating-point-and-non-point-sources-basins).
HSPF is a semi-distributed model that simulates water and contaminant transport over spatially distributed areas within a given catchment [61]. It is the combination of several previous models [62][63][64][65][66] and it is able to simulate continuous dynamic events, of both hydrologic and water quality processes, through surface and groundwater regimes in the catchment. The application of HSPF to model water processes and contaminant transport has been reported in many previous studies and for a wide range of types of catchments and environmental conditions [7,16,17,19,21,67,68].
The calibration process of water quantity and quality was carried out by coupling HSPF with GIS in the delineation of the river catchment data. The calibration and validation parameter values were estimated until good agreement was achieved between observed and simulated flowrates and water quality, which was then compared to the Portuguese legislation for surface freshwater.
Daily and monthly flowrates and quality constituents were calibrated and validated at the hydrometric station "Fragas da Torre", hereafter FT (Figure 1b), for the periods of 1950-1974 and 1980-2011, respectively. This is the only hydrometric station with enough observed data in the catchment to perform adequate model calibration and validation. The daily flowrates were obtained from the Sistema Nacional de Informação de Recursos Hídricos (SNIRH, www.snirh.pt). The verification of the model performance was assessed by three statistical metrics (performance measures): Coefficient of Determination (R 2 ), Percent Bias (PBIAS) and Nash-Sutcliffe efficiency coefficient (E) [69]. The Nash-Sutcliffe model efficiency coefficient can range from −∞ to 1, where 1 is a perfect fit of modelled discharge to the observed data, whereas zero indicates that the model prediction is as accurate as the observed mean.
Regarding the water quality constituents, owing to the lack of timeliness data and accuracy, as well as the low number of records, a different approach was undertaken for calibration. It is not reasonable to expect the model to simulate a daily average concentration equal to an instant observed value, as the recorded values are typically sampled infrequently (monthly at best) and only represent an instant in time. Therefore, observed long-term means were compared against simulated long-term means to reach an agreement for calibration (best-fit parameter set). The monthly variation is also shown as the final simulation result. The calibrated parameters were water temperature, biochemical oxygen demand, dissolved oxygen, nitrates, orthophosphates, chlorophyll-a, faecal coliforms, and total suspended solids. In HSPF, nitrogen and phosphorus are modelled based on nutrient cycle empirical relationships, under the modules of PQUAL (quality constituents for pervious land, PERLND), IQUAL (quality constituents for impervious land, IMPLND), and QGUAL (quality constituents for reaches, RCHRES). Phosphorus load is usually correlated to the sediment load while using a power factor approach, while nitrogen, which is not normally deposited into the sediments, uses a method of accumulation and depletion based on the first-order wash off rate. Biological, chemical, and physical processes in streams were also simulated.

Environmental Variables, Biological and Ecological Indicators
Land-use/Land Cover (LU/LC) data (Table 1) were obtained by GIS from COS 2018 according to DGT (2019) (see Table S2, Supplementary Material). A proximity method was used to overcome the absence of strong associations between environmental descriptors at higher spatial scales and local aquatic communities, defining a circle with a radius of 1 km around each site based on Cortes et al. [70] ( Figure 2). This buffer was imposed on LU/LC maps in order to obtain the percentages of each LU for each site. These variables were calculated from catchment to the river line. Temperature, flow, volume, and river size variables were derived from the above described models (see Table S2, Supplementary Material).

Environmental Variables, Biological and Ecological Indicators
Land-use/Land Cover (LU/LC) data (Table 1) were obtained by GIS from COS 2018 according to DGT (2019) (see Table S2, Supplementary Material). A proximity method was used to overcome the absence of strong associations between environmental descriptors at higher spatial scales and local aquatic communities, defining a circle with a radius of 1 km around each site based on Cortes et al. [70] (Figure 2). This buffer was imposed on LU/LC maps in order to obtain the percentages of each LU for each site. These variables were calculated from catchment to the river line. Temperature, flow, volume, and river size variables were derived from the above described models (see Table S2, Supplementary Material).  Physicochemical in situ parameters were measured in 20 sampling sites, which were distributed along the Paiva River catchment ( Figure 2). Water samples were also collected, preserved in cooling conditions (4 • C), and then transported to the laboratory for analyses using standard methods. The analysed parameters were alkalinity (mg HCO 3 , and chromium (µg Cr·L −1 ). The values were classified according to the criteria of the ex-Portuguese Water Institute (INAG). This classification comprises five different classes and a site is considered to be unimpaired if it belongs to class A and impaired if it belongs to classes D or E. The environmental variables were standardized prior to all analyses in order to preserve the original scale. Invertebrate assemblages were simultaneously collected with physicochemical parameters for each site with a hand-net following a WFD compliant protocol (INAG I.P., 2008). The samples were sorted and preserved with an ethanol concentration of ca. 70%. Subsequently, they were keyed to the family level, except for individuals belonging to the Oligochaeta and Hydracarina classes. The work was carried out during the autumn of 2018.
The macroinvertebrate data were organized into two data groups: (a) relative abundance of taxa, resulting in a matrix comprising 52 invertebrate families (see Table S3, Supplementary Material). The data of taxa were log(x+1) transformed before subsequent analyses; (b) biological traits (food and feeding habits), ecological traits (longitudinal and transversal distribution-zonation, current and microhabitat preferences, locomotion type), and conventional metrics (richness, composition and dominance measures, descriptors, and biotic indices) reflecting the life history of taxa and adaptive responses to environmental factors (see Table S4 [72], Supplementary Material). To obtain this second matrix, macroinvertebrate data were imported into the ASTERICS 4.03 software (AQUEM/STAR Ecological River Classification System) in order to calculate traits and metrics. The data were log(x + 1) transformed, except for metrics expressed in percentage. These metrics were not transformed, since they were normally distributed (see Table S4, Supplementary Material).
The Northern Portugal Invertebrate Index (IPtI N ) was used in order to assess the biological status of freshwater [73][74][75][76][77]. This index depicts information on how to differentiate the communities that live in the river sections with good ecological integrity. Its calculation is based on the presence of benthic invertebrates in freshwater due to their sensitivity to a myriad of contamination sources [1,78]. The ecological status was determined with the software Amiib@-"Application for the calculation of Benthic Invertebrate Metrics and Indices" developed by the Portuguese Water Institute (Lisbon, Portugal). It allows for the determination of several Portuguese official indexes for assessing the ecological status of a river, based on the benthic invertebrate biological element, including IPtI N .

Correspondence Analysis
Ordinations of sampling sites according to macroinvertebrate relative abundance and traits/metrics were configured in order to explore the ecological patterns (i.e., whether or not the axes of ordinations could characterize gradients) by using correspondence analysis (CA), starting with 25 environmental variables. CAs were accomplished for both biological datasets, with the environmental matrices separately (physicochemical and LU/LC) providing bi-plots of all datasets. This multivariate analysis was applied in order to reduce the dimensionality and extract the dominant factors. CA can be viewed as a method projecting data into a multidimensional space, where the chi-square distance metric is preserved, but also as a method that computes latent predictor variables for (approximated) models of unimodal response. The CA was carried out using CANOCO 5 (version 5.12, Microcomputer Power, Ithaca, NY, USA) [79].
Relating species traits/metrics to physicochemical and LU/LC can provide important insights into the structure and functioning of stream communities. The Spearman correlation test was used to avoid the collinearity within traits/metrics. Hence, a non-redundant matrix of 13 traits/metrics from an initial matrix with 77 was built. The remaining traits belong to current preferences, food and feeding habits, locomotion type, and taxonomic groups (number of taxa) (Table S2, Electronic Supplementary Material). The main objectives are to demonstrate the importance of the different macroinvertebrate traits that are related to perturbation patterns and identify which of them allow more effective and accurate use in future monitoring works.

Current Conditions
The mean annual temperature for the period 1950-2015 ( Figure 3a) shows higher temperatures along the river line, reaching values as high as 15.5 • C, while much lower temperatures can be found in the mountains north and south of the river (~11 • C), thus showing a clear footprint of elevation (negative vertical temperature gradient). For mean temperature (see Figure S1, Supplementary Material), the observed variation is between 20 • C in summer and 5.5 • C in winter. Hence, according to the Köppen-Geiger climate classification, the Paiva River catchment presents a Csb climate (Mediterranean climate with cool summers), or predominant mesomediterranean climate with supramediterranean climate at the highest elevations [80,81]. Regarding the mean annual precipitation for the period 1950-2015 (Figure 3d), it increases as we move upstream the river (from 1050 to 1200 mm), but it shows higher values in the mountains (1400 mm). The precipitation pattern is also a clear footprint of the orography, as the condensation barrier effect that is imposed by the mountains to the Atlantic westerly and moist winds underlies the observed precipitation pattern, which is further enhanced by the Atlantic-facing deep river canyon, with a tunnelling and uplift effect on the air masses. During winter (DJF), the mean precipitation in the Paiva River catchment shows values of 620 mm, whereas much lower values (ca. 85 mm) are recorded in summer (JJA) (see Figure S1, Supplementary Material). In terms of aridity, the climate can be considered to be humid everywhere, as the annual mean precipitation is higher than annual mean potential evapotranspiration [82]. Mid-seasons (MAM and SON) show similar spatial patterns, with a variation of temperature between 10.5 and 15.5 • C and precipitation between 275 and 425 mm.
Although the Paiva River catchment depicts a typical Mediterranean climate (see Figure S2, Supplementary Material), the maritime influence from the North Atlantic is apparent, with a commonly short dry season (mid-June to mid-August) and considerably wet winters. While, in winter, frequent synoptic-scale disturbances travel from the North Atlantic to the region, thus generating high precipitation amounts, the prevalence of strong and northerly displaced Azores high-pressure systems in summer, accompanied by air sinking, is largely unfavourable to rain-generating conditions and it explains the summer minimum [83]. Therefore, the strong seasonality in the precipitation regime has significant effects on river flow, as well as on water resources and water quality. The irregularity in the monthly precipitation amounts is also noteworthy (see Figure S2, Supplementary Material).

Climate Change Projections
The projections for both RCP scenarios (for the period 2041-2070) (Figure 3b,c) show an overall increase of temperature for the whole of the region, where scenario RCP8.5 shows an increase for the entire catchment up to 2 • C. Concerning the precipitation scenarios for the period 2041-2070 (Figure 3e,f), there is a projected decreased in precipitation for scenario RCP4.5 of approximately 40 mm (maximum), while, in RCP8.5, there is an observed increase of precipitation in the eastern part of the catchment, varying from −40 to 30 mm. 620 mm, whereas much lower values (ca. 85 mm) are recorded in summer (JJA) (see Figure S1, Supplementary Material). In terms of aridity, the climate can be considered to be humid everywhere, as the annual mean precipitation is higher than annual mean potential evapotranspiration [82]. Midseasons (MAM and SON) show similar spatial patterns, with a variation of temperature between 10.5 and 15.5 °C and precipitation between 275 and 425 mm. Although the Paiva River catchment depicts a typical Mediterranean climate (see Figure S2, Supplementary Material), the maritime influence from the North Atlantic is apparent, with a commonly short dry season (mid-June to mid-August) and considerably wet winters. While, in winter, frequent synoptic-scale disturbances travel from the North Atlantic to the region, thus The chronograms of the area-mean climate variables (Figure 4), averaged over the whole catchment, show an increasing trend for temperature and a decreasing for precipitation for both RCP 4.5 and RCP 8.5 scenarios in the Paiva catchment. The projected temperature shows an increase of 0.5 • C and 1.1 • C under RCP 4.5 and RCP 8.5, respectively. While projected precipitation decreases when compared to the historical period. Overall, there is an average annual decrease of 230 mm and 214 mm of precipitation for scenarios RCP 4.5 and RCP8.5, respectively.
Two indices of precipitation extremes were also computed in order to better understand the precipitation regime of the Paiva River catchment and its likely changes in the future. The number of consecutive dry days (CDD) and the relative contribution of precipitation above the 95-percentile to total precipitation (R95PTot) ( Figure 5) are showed here as indicators to assess seasonal droughts and extreme precipitation events. Both of the indicators were calculated for the entire period of the gridded dataset  and the period 2041-2070 for both RCP scenarios. CDD shows a similar spatial pattern as the seasonal precipitation, with 105 consecutive dry days upstream, to 85 days, downstream. On the contrary, R95PTot shows that the severest precipitation events tend to occur at higher altitudes (inner part of the catchment). In future scenarios, CDD tends to decrease and the extremes of precipitation tend to be more frequent and/or intense, particularly in the upper area of the catchment.

Climate Change Projections
The projections for both RCP scenarios (for the period 2041-2070) (Figure 3b,c) show an overall increase of temperature for the whole of the region, where scenario RCP8.5 shows an increase for the entire catchment up to 2 °C. Concerning the precipitation scenarios for the period 2041-2070 ( Figure  3e,f), there is a projected decreased in precipitation for scenario RCP4.5 of approximately 40 mm (maximum), while, in RCP8.5, there is an observed increase of precipitation in the eastern part of the catchment, varying from −40 to 30 mm.
The chronograms of the area-mean climate variables (Figure 4), averaged over the whole catchment, show an increasing trend for temperature and a decreasing for precipitation for both RCP 4.5 and RCP 8.5 scenarios in the Paiva catchment. The projected temperature shows an increase of 0.5 °C and 1.1 °C under RCP 4.5 and RCP 8.5, respectively. While projected precipitation decreases when compared to the historical period. Overall, there is an average annual decrease of 230 mm and 214 mm of precipitation for scenarios RCP 4.5 and RCP8.5, respectively. Two indices of precipitation extremes were also computed in order to better understand the precipitation regime of the Paiva River catchment and its likely changes in the future. The number of consecutive dry days (CDD) and the relative contribution of precipitation above the 95-percentile to total precipitation (R95PTot) ( Figure 5) are showed here as indicators to assess seasonal droughts and extreme precipitation events. Both of the indicators were calculated for the entire period of the gridded dataset  and the period 2041-2070 for both RCP scenarios. CDD shows a similar spatial pattern as the seasonal precipitation, with 105 consecutive dry days upstream, to 85 days, downstream. On the contrary, R95PTot shows that the severest precipitation events tend to occur at higher altitudes (inner part of the catchment). In future scenarios, CDD tends to decrease and the extremes of precipitation tend to be more frequent and/or intense, particularly in the upper area of the catchment.

Flowrates
Based on the records from the FT station, the medians of the daily mean flowrates vary from 25.5 m 3 s −1 , in January, to 1.16 m 3 s −1 , in August, while daily mean flowrates range from 0.19 to 920 m 3 s −1 . Hence, the observed daily flowrates hint at strong seasonality, peaking in winter, with very pronounced maximum flowrates, which are a typical feature of the Mediterranean-type flow regimes and they reflect the aforementioned precipitation regime. The values of the hydrologic parameters for the model calibration are within the range of those presented in the literature [84].
The HSPF model performed well with respect to the simulation of river flow in the Paiva River (Figure 6), since the fits between observed and simulated flows present high Nash-Sutcliffe coefficients (E) values, as well as high coefficients of determination (R 2 ), for both calibration and validation, as shown hereafter. The deviation of volumes (Dv) also confirms the goodness-of-fit of the model, with a deviation of 3% and −11% for calibration and validation, respectively, though, for validation, the model slightly underestimates the flowrates. The Nash-Sutcliffe model efficiency coefficient shows a value of 0.95 for calibration and 0.76 validation. Based on previously published recommendations [85], model simulations with values of E above 0.75 are generally considered to be satisfactory. The coefficients of determination also confirm the goodness-of-fit of the model, with particularly good performance >0.95 for calibration and 0.82 for validation. It is important to consider many factors when modelling water quality in a catchment, of which the quality of the observed data is of the utmost importance. In the Paiva River catchment, there is considerable variability of the observed values and a general lack of homogenous and continuous data, which is an important limitation to the analysis.
Concerning the distribution of the projected precipitation on the daily timescale (Figure 7), an increase in skewness under both RCP are found (more extremes), despite the aforementioned decrease in the annual means. These changes are in general agreement with Figure 5. These changes combined lead to a decrease of flowrate variability (distributional dispersion) from the historical period to RCP4.5 (lower precipitation), and an increase of the flowrate variability from scenario RCP4.5 to RCP8.5 (more extremes).

Flowrates
Based on the records from the FT station, the medians of the daily mean flowrates vary from 25.5 m 3 s −1 , in January, to 1.16 m 3 s −1 , in August, while daily mean flowrates range from 0.19 to 920 m 3 s −1 . Hence, the observed daily flowrates hint at strong seasonality, peaking in winter, with very pronounced maximum flowrates, which are a typical feature of the Mediterranean-type flow regimes

Water Quality
In this study, the water quality constituents were simulated on a daily timescale by HSPF, while observed data were only available at the monthly timescale and with numerous gaps over the simulation period. Therefore, the calibration and validation of the model were performed in order to match the observed long-term mean against the simulated long-term mean, while displaying the monthly simulated hydrographs (Figure 8) to understand the temporal evolution and seasonality of the quality constituents.
validation, as shown hereafter. The deviation of volumes (Dv) also confirms the goodness-of-fit of the model, with a deviation of 3% and −11% for calibration and validation, respectively, though, for validation, the model slightly underestimates the flowrates. The Nash-Sutcliffe model efficiency coefficient shows a value of 0.95 for calibration and 0.76 validation. Based on previously published recommendations [85], model simulations with values of E above 0.75 are generally considered to be satisfactory. The coefficients of determination also confirm the goodness-of-fit of the model, with particularly good performance >0.95 for calibration and 0.82 for validation. It is important to consider many factors when modelling water quality in a catchment, of which the quality of the observed data is of the utmost importance. In the Paiva River catchment, there is considerable variability of the observed values and a general lack of homogenous and continuous data, which is an important limitation to the analysis. Concerning the distribution of the projected precipitation on the daily timescale (Figure 7), an increase in skewness under both RCP are found (more extremes), despite the aforementioned decrease in the annual means. These changes are in general agreement with Figure 5. These changes combined lead to a decrease of flowrate variability (distributional dispersion) from the historical

Water Quality
In this study, the water quality constituents were simulated on a daily timescale by HSPF, while observed data were only available at the monthly timescale and with numerous gaps over the simulation period. Therefore, the calibration and validation of the model were performed in order to Several parameters were adjusted first for calibration based on the literature [7,84,[86][87][88]. The monthly nutrient concentration in interflow (MON-IFLW-CONC, mg·L −1 ) and groundwater (MON-GRND-CONC, mg·L −1 ) had significant impacts on phosphorus and nitrogen yields, followed by the ratio of constituent yield to sediment outflow (POTFS), the rate of surface runoff that will remove 90% of the quality constituent (WSQOP), and the maximum storage of the quality constituent that can be stored on the land surface (SQOLIM). For nitrates (NO 3 ) calibration, the observed and simulated concentrations were 2.257 and 2.218 mg·L −1 , respectively. For validation, those values were 2.200 and 1.987 mg·L −1 . These results indicate a good agreement between the observed and simulated NO 3 concentrations, with differences of 1.7% for calibration and −9.7% for validation. There is an underestimation of NO 3 concentrations in the validation period, although most of the records of the observed values from the water quality station stated less than 2. Moreover, the concentration of NO 3 for the entire period is below the maximum recommended value (MRV) (25 mg·L −1 ) of surface freshwater intended for the production of drinking water for type A1 treatment (physical treatment and disinfection) [89]. Regarding the orthophosphates (PO 4 3  (Table 2), except for water temperature (WT) and dissolved oxygen (DO). The remaining parameters show simulated values that are below the established MRV, except for BOD in the calibration period, depicting three occurrences that overlap the limit (3 mg O 2 L −1 ), and for faecal coliforms (FC) across the entire period, with values that are above 20 UFC 100 mL −1 . The faecal coliform concentrations are below 2000 UFC mL −1 , meaning that the type A2 treatment (physical and chemical treatment plus disinfection) would be necessary to comply with the above legislation. The percentage difference values (PBIAS) for calibration and validation of the water quality parameters show good performance according to the literature [85]. Furthermore, a good calibration of quality constituents is not only dependent on the frequency, timing, and the total number of samples, but also of good hydrology calibration. In the particular case of the Paiva River, even though a good hydrology calibration was achieved (based on the criteria results) at FT, it is more difficult to reach a meaningful water quality calibration, since there is only one hydrometric station with enough data for model calibration. This may explain why the HSPF results for validation showed a consistent bias towards an underestimation of the concentrations (PBIAS < 0) of most of the quality constituents, as the validation of the hydrological model also showed a deviation between the simulated and observed flowrate of −11%. Further, the observed quality constituent records are infrequent and/or below the detection limit of the legislated methodology [89].

Current Ecological Status versus Water Quality
The macroinvertebrate sampling revealed a total of 46262 specimens within 87 families. The most representative families were: Simuliidae (9860-21.7%, present in all sampling sites), Chironomidae (5688-12.5%), Philopotamidae (5449-12.0%), Elmidae (4378-9.6%), Baetidae (3550-7.8%), Hydropsychidae (3039-6.7%), Heptageniidae (2035-4.5%), and Caenidae (1034-2.3%) (Table S3). These eight families represent 77.1% of the all captured individuals and 43.2% of them pertain to the orders Ephemoptera (mayflies) and Trichoptera (caddisflies) that are pollution sensitive [90][91][92]. Further, Simuliidae is considered to be a reliable indicator of good water quality [93,94]. Notwithstanding, Baetidae, Caenidae (Ephemeroptera), and Hydropsychidae (Trichoptera) families exhibited tolerance to organic pollution and habitat disturbance [95,96]. According to Raburu et al. [97] and Masese et al. [98], this tolerance enables them to be broadly distributed, with higher abundance in areas with this kind of pollution. Despite this fact, it can be highlighted that the majority of sampling sites (65%) have more than 40% of EPT-Taxa [%] (min-max: 4.5%-79.0%) and all of them have more than three very intolerant EPT families (Table S3). Additionally, in agreement with Chirhart [99], the presence of moderate numbers of intolerant taxa is an indicator of good aquatic health. This corroborates the general good biological status that was verified throughout the Paiva catchment ( Figure 9) using the IPtI N index. According to Marshall et al. [100], this same author, biotic indices around the family level of macroinvertebrates preclude the need for species-level identification, which can be time consuming and requires taxonomically trained specialists.

Current Ecological Status Versus Water Quality
The macroinvertebrate sampling revealed a total of 46262 specimens within 87 families. The most representative families were: Simuliidae (9860-21.7%, present in all sampling sites), Chironomidae (5688-12.5%), Philopotamidae (5449-12.0%), Elmidae (4378-9.6%), Baetidae (3550-7.8%), Hydropsychidae (3039-6.7%), Heptageniidae (2035-4.5%), and Caenidae (1034-2.3%) (Table S3). These eight families represent 77.1% of the all captured individuals and 43.2% of them pertain to the orders Ephemoptera (mayflies) and Trichoptera (caddisflies) that are pollution sensitive [90][91][92]. Further, Simuliidae is considered to be a reliable indicator of good water quality [93,94]. Notwithstanding, Baetidae, Caenidae (Ephemeroptera), and Hydropsychidae (Trichoptera) families exhibited tolerance to organic pollution and habitat disturbance [95,96]. According to Raburu et al. [97] and Masese et al. [98], this tolerance enables them to be broadly distributed, with higher abundance in areas with this kind of pollution. Despite this fact, it can be highlighted that the majority of sampling sites (65%) have more than 40% of EPT-Taxa [%] (min-max: 4.5%-79.0%) and all of them have more than three very intolerant EPT families (Table S3). Additionally, in agreement with Chirhart [99], the presence of moderate numbers of intolerant taxa is an indicator of good aquatic health. This corroborates the general good biological status that was verified throughout the Paiva catchment ( Figure 9) using the IPtIN index. According to Marshall et al. [100], this same author, biotic indices around the family level of macroinvertebrates preclude the need for species-level identification, which can be time consuming and requires taxonomically trained specialists. Overall, the quality is excellent, apart from the sampling sites number 4, 16, 18, and 20 ( Figure  9), which reveal good ecological status. This relatively good quality had also been reported by Sousa et al. (2013) [76] when evaluating the ecological status of a Margaritifera margaritifera population at the Southern Edge of its Distribution, i.e., in the Paiva River.
We performed CA analyses to test the likely influence of environmental variables at different scales (LU/LC and water quality+flow) on the distribution of macroinvertebrate species and traits/metrics ( Figure 10). Overall, the quality is excellent, apart from the sampling sites number 4, 16, 18, and 20 (Figure 9), which reveal good ecological status. This relatively good quality had also been reported by Sousa et al. (2013) [76] when evaluating the ecological status of a Margaritifera margaritifera population at the Southern Edge of its Distribution, i.e., in the Paiva River.
We performed CA analyses to test the likely influence of environmental variables at different scales (LU/LC and water quality+flow) on the distribution of macroinvertebrate species and traits/metrics ( Figure 10).
The invertebrate taxonomic data had a gradient of 1.4 SD units along the main axis of the CA (Figure 10). Seven variables were selected from the initial number of 25 environmental variables at reach scale, following a forward selection approach (p < 0.05). The environmental variables, which included water quality+flow, explained 45.92% of the variance (with a total variation of 0.67967). The first two axes of the CA explained 22.43% and 16.05% of the total variance, respectively.
CAs that were accomplished in the species matrix and water quality+flow variables revealed in Axis-2 a very clear impact gradient related to chemical and organic contamination (Figure 10a).
On the top left of the ordination space are sites (P2, P7, P14) impacted by organic contamination related with the increase of nitrates. This situation is a consequence of intensive eucalyptus afforestation observed mainly in the western part of river catchment (P2 and P7). Conversely, on the bottom, sites P1, P4, P5, P8 and P19 are associated with chemical contamination namely high concentrations of Al, Ni, and Zn, although the values obtained respect the minimum environmental quality objectives according Portuguese legislation. According to Rosseland et al. [101], aluminium acts as a toxic agent on gill-breathing animals such as fish and invertebrates, by causing loss of plasma-and haemolymph ions leading to osmoregulatory failure. Yet, Sparling et al. [102] refer that in general, aquatic invertebrates are less sensitive to Al toxicity and acidity than fish. Regarding zinc, Carvalho et al. [103] report that, even in relatively low concentrations, zinc metals can alter the reproductive pattern of freshwater invertebrates, compromising the aquatic ecosystem and its recovery capacity. Relatively to Ni, Eisler [104] state that its concentrations in unimpaired aquatic ecosystems are typically below 10 µg/L. CA analyses show that 50% of the studied sites (P1, P3, P4, P5, P6, P15, P16, P17, P18, P19) had a relationship with this heavy metal, though only detected in P4, P16, and P18 where good biological quality was achieved. While Ni is considered to be essential for a wide variety of animals' species [105], several studies revealed a Ni-related depression of immune system in both vertebrates and invertebrates [104,106]. According the Portuguese Mineral Resources and Occurrence Information System (SIORMINP), provided by National Energy and Geology Laboratory, I.P, these three metallic minerals occur in the studied area. Although the ecological status of the catchment is very good to excellent, according to the IPtI N index, CA analyses revealed minor contamination and established relationships between nutrients and macroinvertebrates.
CAs that were made with species matrix and LU/LC variables resulted in the selection of five LU/LC variables from the initial number of 11, following a forward selection approach (p < 0.05). This analysis showed a gradient related to agriculture in the Axis-1 (Figure 10b). The LU/LC variables only explained 24.59% of the variance (with a total variation of 0.67967). Ordination of sites on the CA plot can be interpreted by the presence of species in each assemblage. Larval Ephemerellidae (mayfly), Gomphidae (dragonflies), Chloroperlidae and Perlidae (stoneflies), Dytiscidae and Scirtidae (beetles), and Corbiculidae (clam) were the most influential taxa on Axis-1 (Figure 10b). Although this gradient that opposes agriculture and shrubs to the presence of eucalyptus (P1-P5 and P7) and urban areas (P13 and P18) has been verified, to notice, however, that both situations present intolerant families to pollution, reflecting the good biological quality obtained in the studied area. This may be justified by agricultural activities, mostly of extensive nature, and with the presence of small urban areas. In what concerns eucalyptus plantations [107], still the consequences of these plantations for the functioning of water courses are largely unknown. Fundazioa [107] refers that, despite the decomposition of its dead leaves having fewer problems in the initial time, due to the phenolic components, this material tends to be used in rivers as much as that from deciduous trees like autochthonous species. Additionally, Abelho and Graça [108], when comparing the dynamic of the organic matter (eucalyptus versus chestnut leaf litter), available for the aquatic macroinvertebrates, found that the monoculture of eucalyptus affects low order stream communities. Notwithstanding, it was verified that the impact may be attenuated if natural riparian vegetation is kept in those areas of plantation.
Finally, Figure 10d shows summary results of taxa diversity within samples. The CA highlight that the higher biodiversity was recorded in River Paiva tributaries (P2, P7, P11, and P14), which present excellent quality, linked with well-preserved areas.
The traits/metrics data had a gradient of 0.8 SD units along the main axis of the CA (Figure 11). From the initial number of 25 environmental variables at reach scale, following a forward selection approach (p < 0.05), six variables were selected. Water quality+flow, explained 35.54 % of the variance (with a total variation of 0.15795). The first two axes of the CA explained 30.45% and 23.15% of the total variance, respectively. Contrary to what was verified with the analyses of invertebrate taxonomic data, the traits/metrics gradient organic versus chemical disturbance was reflected in Axis-1 (Figure 11a). On the bottom left of the ordination space are sites (P2, P4, P7, P11, P12, P14), impacted by organic contamination also related with the increase of nitrates. On the right site of the ordination, sites P6, P13, P16, P18, P19, and P20 are associated with chemical contamination, namely high concentrations of Mn and Zn. Manganese also respects the minimum environmental quality objectives according Portuguese legislation. Manganese, which is an essential element actively assimilated and utilized by both plants and animals, can be significantly bio concentrated by aquatic biota at lower trophic levels [109]. According to this same author, freshwater mollusks and crustaceans are the most manganese-sensitive freshwater invertebrates, followed by arthropods and oligochaetes. Oligochaetes and gastropods are in the left part of the Axis-1 corroborating this idea, as it can be seen in the Figure 11a. Heal [110] refers that afforestation of upland areas increase manganese concentrations in surface waters. Additionally, has been widely documented that catchments planted with conifers that acidify soil and water are associated with enhanced manganese concentrations in surface waters [110][111][112]. Thus, by observing Figure 2 the coniferous woodland in mainly located at the upper part of Paiva catchment which is in agreement with the previous statement. As mentioned before, even though the ecological status is good to excellent, the CA analyses for the traits/metrics also demonstrate a correlation between species and organic contamination.  (Figure 10). Seven variables were selected from the initial number of 25 environmental variables at reach scale, following a forward selection approach (p < 0.05). The environmental variables, which included water quality+flow, explained 45.92% of the variance (with a total variation of 0.67967). The . Only significant and independent environmental variables and more fitted taxa are shown. The red arrow shows the direction of the disturbance gradient and the green arrow the biodiversity gradient, common to all diagrams. Sites, taxa, and environmental codes are given in Table 1 and Table S2-Electronic Supplementary Material. catchments planted with conifers that acidify soil and water are associated with enhanced manganese concentrations in surface waters [110][111][112]. Thus, by observing Figure 2 the coniferous woodland in mainly located at the upper part of Paiva catchment which is in agreement with the previous statement. As mentioned before, even though the ecological status is good to excellent, the CA analyses for the traits/metrics also demonstrate a correlation between species and organic contamination. The ordination made with LU/LC variables resulted in the selection of five (Eu, UR, CW, Ag, and SH) variables from the initial number of 11, following a forward selection approach (p < 0.05). A gradient related to agriculture in the Axis-1 (Figure 11b). This analysis showed a separation between the sites associated with the areas of large eucalyptus plantations (monoculture) against the remaining LU/LC in which agriculture has the greatest influence. The LU/LC variables only explained 25.03% of the variance (with a total variation of 0.15795). Shredders life-history characteristics are linked with agriculture and shrubs areas that mainly occur in the upper part (P19 and P20) and lower western part (P6) of River Paiva catchment. Masese et al. [98] stated that shredders are clearly very important in the transformation process of organic matter in forested and deforested streams. In our work, despite the shredders being associated with the presence of (c) CA plot with the 20 sampling sites from the Paiva catchment. Only significant and independent environmental variables and more fitted traits are shown. The red arrows show the direction of the disturbance gradient and are common to all diagrams. Sites, taxa and environmental codes are given in Table 1 and Table S2-Electronic Supplementary Material.
The ordination made with LU/LC variables resulted in the selection of five (Eu, UR, CW, Ag, and SH) variables from the initial number of 11, following a forward selection approach (p < 0.05). A gradient related to agriculture in the Axis-1 (Figure 11b). This analysis showed a separation between the sites associated with the areas of large eucalyptus plantations (monoculture) against the remaining LU/LC in which agriculture has the greatest influence. The LU/LC variables only explained 25.03% of the variance (with a total variation of 0.15795). Shredders life-history characteristics are linked with agriculture and shrubs areas that mainly occur in the upper part (P19 and P20) and lower western part (P6) of River Paiva catchment. Masese et al. [98] stated that shredders are clearly very important in the transformation process of organic matter in forested and deforested streams. In our work, despite the shredders being associated with the presence of agriculture, contrary to the huge amount of studies that prove the decrease of these organisms with this type of uses [113][114][115][116][117], it should be noted that these areas are also associated with the strong presence of scrub and broadleaf woodlands. These areas contribute with great amount of leaves, their main source of food, for breakdown processes.
As verified by Larrañaga et al. [118], eucalypt plantations in a northern Spain stream has a negative effect on macroinvertebrate communities mostly on the density of shredders. This is in accordance with what was observed in Paiva catchment, where shredders were put in the opposite end in the ordination space (Axis-1) to the monoculture of eucalyptus. The burrowing/boring traits (P2, P3, and P4) were associated with this LU/LC. The presence of invertebrates that burrows in the substrate could be favored by the increase of the sedimentation in river channels [119]. Thus, it is precisely in this part, further downstream of the Paiva basin, where a greater amount of finer sediment is found in the riverbed as a consequence of the afforestation and/or harvesting of the eucalyptus monoculture that requires intervention in the soil with relatively short periods due to the exploitation of timber.
Swimming/skating were dominant in urban and coniferous areas (P13, P16, and P18). The potential for land use effects on the composition of aquatic invertebrate communities to alter the composition of behavioral related traits was also seen in research by Wang et al. [116], which found that forested streams exhibited higher relative abundances of Swimmer taxa when compared with the agricultural stream.

Conclusions
The HSPF model was applied to a 796 km 2 catchment in northern Portugal, the Paiva River catchment, a predominantly forest and agricultural catchment, aiming to validate its performance in simulating hydrologic and water quality responses. The flowrate and water quality components were calibrated against the observed data that were collected at the hydrological station of "Fragas da Torre" and extrapolated for the entire watershed. Overall, the runoff volumes predicted by HSPF agreed well with the observed data at the monitoring site within the catchment. Both model calibration and validation show satisfactory criteria values. Furthermore, because simulated water quality processes are based on hydrologic processes, the quality of the hydrology calibration has a major impact on the quality of the calibration of the constituent concentrations. Nonetheless, a satisfactory calibration was achieved. The status of water quality is considered to be good when compared with the maximum recommended value for surface freshwater intended for the production of drinking water. According to the HSPF model output results, only BOD, PO 4 , and faecal coliforms would require a secondary treatment to comply with the legislation in force. Regarding the ecological status, an overall excellent status is observed in the Paiva River catchment, with only four sites disclosing good ecological status.
Climate change is likely to have extensive effects on water quality, in the particular case of Paiva River catchment, due to increase of extreme events in both precipitation and temperature, both of which will affect the concentration of water pollutants and erosion, especially for scenario RCP8.5. The simulated flowrates for future scenarios show an increasing trend that may contribute to a better water quality status (providing the current anthropogenic pressures), thus reflecting a continuous likelihood of good ecological status.
Correspondence analyses highlight some potential predictors of macroinvertebrate populations, which can be used to implement population distribution models in the Paiva River catchment. These models can be of foremost relevance for developing projections for macroinvertebrates under future regional climate change scenarios.
While the present study outcomes refer to a specific watershed, they may be generalized to other watersheds in Portugal and in the Mediterranean basin, as the climate change projections for these regions indicate future drier and warmer climates.
Nature-based solutions can offer sustainable approaches in order to cope with climate change mitigation and adaptation challenges. Namely, nature-based climate adaptation aims to preserve ecosystem services and reduce the impact of the projected negative effects of climate change (e.g., floods, droughts, and heat waves).  Table S1: Integrated General Circulation Models (GCMs) with Regional Climate Models (RCMs) used to assess climate change impacts on Paiva River hydrology, Table S2: Environmental variables (water quality and LU/LC), Code and units used in the analysis, Table S3: Invertebrate taxa, respective codes and expected response to disturbance (based in Alba Tercedor, et al. and Cabecinha, et al., Table S4: Invertebrate Traits, respective codes, transformation and expected response to disturbance (based in OLIVEIRA; https://www.freshwaterecology.info/abbreviations.php). In Bold are represented the selected traits/metrics that were used in the Correspondence analyses (CA).

Author Contributions:
In the overall context of this paper, the contributions are as follow: climate assessment done by A.R.F. and J.A.S., Hydrological model by A.R.F., macroinvertebrates sampling and water analyses done by A.R.F., E.C., S.G.P.V. and S.M.M., data analyses by E.C., J.L.M., R.M.V.C. and S.G.P.V. All authors provided critical feedback and helped shape the research, analysis and manuscript. All authors have read and agreed to the published version of the manuscript.