Assessment of Ecological and Hydro-Geomorphological Alterations under Climate Change Using SWAT and IAHRIS in the Eo River in Northern Spain

Magnitude and temporal variability of streamflow is essential for natural biodiversity and the stability of aquatic environments. In this study, a comparative analysis between historical data (1971–2013) and future climate change scenarios (2010–2039, 2040–2069 and 2070–2099) of the hydrological regime in the Eo river, in the north of Spain, is carried out in order to assess the ecological and hydro-geomorphological risks over the short-, medium- and long-term. The Soil and Water Assessment Tool (SWAT) model was applied on a daily basis to assess climate-induced hydrological changes in the river under five general circulation models and two representative concentration pathways. Statistical results, both in calibration (Nash-Sutcliffe efficiency coefficient (NSE): 0.73, percent bias (PBIAS): 3.52, R2: 0.74) and validation (NSE: 0.62, PBIAS: 6.62, R2: 0.65), are indicative of the SWAT model’s good performance. The ten climate scenarios pointed out a reduction in rainfall (up to −22%) and an increase in temperatures, both maximum (from +1 to +7 °C) and minimum ones (from +1 to +4 °C). Predicted flow rates resulted in an incrementally greater decrease the longer the term is, varying between −5% (in short-term) and −53% (in long-term). The free software IAHRIS (Indicators of Hydrologic Alteration in Rivers) determined that alteration for usual values remains between excellent and good status and from good to moderate in drought values, but flood values showed a deficient regime in most scenarios, which implies an instability of river morphology, a progressive reduction in the section of the river and an advance of aging of riparian habitat, endangering the renewal of the species.


Introduction
River systems provide ecological services of utter importance both to Earth's biodiversity and society's development [1]. Therefore, adequate management demands a detailed analysis of their structure and behavior. Directive 2000/60/EC of the European Parliament and of the Council establishes a framework for community action in the field of water policy [2]. One of the main goals is to achieve the good status of water bodies, both surface water and groundwater, protecting them and avoiding their deterioration. Their status will be determined by the sum of ecological and chemical characteristics. The assessment of ecological status is evaluated by the combination of biological, physical-chemical and hydro-geomorphological indexes. The latter conditions can be grouped into four categories: physical habitat assessment, riparian habitat assessment, morphological assessment and assessment of hydrological regime alteration [3].
Flow regime has a core role in the river ecosystem, setting up both riparian and aquatic environments and providing habitats with various and dynamic interactions. The full range of intraand inter-annual hydrologic regimes, as well as the magnitude, frequency, timing, duration, and rate of change of water flows, are essential for the natural biodiversity and stability of aquatic environments [4]. Accordingly, the success of conservation in biodiversity and river services depends on our ability to protect and restore the main environmental aspects of the natural flow regime [5][6][7]. Flow regime has a significant influence on the biotic composition, structure and behavior of river ecosystems [8] and performs in all temporal and spatial levels [9]. Bunn and Arthington [7] distinguish three main mechanisms of interactions between hydrological regime and biological communities: flow as the major determinant of habitat, strategies for adapting to natural variability of hydrological regime developed by river species, and patterns of longitudinal and transversal connectivity based on hydrological processes. However, the use of water resources has been carried out, ignoring these processes and self-regulating conditions [10]. The alteration of the hydrologic regime is a key factor that threatens these ecosystems, which, furthermore, interacts with other anthropic disturbances [11]. Natural flow regime may be modified by pressures, such as construction of dams [12], surface [13] and underground water [14] withdrawals and modifications of soil uses [15], whose consequences can be aggravated by the effects of climate change [16].
Indeed, in many regions, due to climate change, it is currently clear that "changing precipitation or melting snow and ice are altering hydrological systems, affecting water resources in terms of quantity and quality" [17]. The Intergovernmental Panel on Climate Change (IPCC) suggested that the Earth is going through a rapid period of climate change and profound changes in climatic variables are expected by the end of the twenty-first century [17][18][19][20]. It is assumed, with some uncertainty, that human actions had a dominant influence on global average temperature change in the 1950-2010 period [21]. In many regions worldwide, hydrologic alterations caused by climate change, such as increased frequency and intensity of floods and droughts, have been observed, causing important human and material damage [22,23]. Numerous studies using combined models of climate change and atmospheric circulation have been carried out [24,25] and they have agreed in predicting a rise in global average temperature and an increase in both global precipitation and evapotranspiration, causing a reduction in the moisture of the soil [26,27] and a significant decrease in runoff [28][29][30]. All of this is compounded by the disruption of habitats and species by mankind and the degradation and loss of the hydrological integrity of river basins [31]. The use of hydrological models plays an essential role, both to identify the current effects of climate change and to assess future scenario results that enable us to take suitable corrective actions. Most studies currently use general circulation model (GCM) simulations to obtain future climatic data. A GCM is a mathematical model of the general circulation of a planetary atmosphere or ocean that attempt to simulate the Earth's climate system. GCMs have unique spatial resolutions, variance magnitudes and processes, so data provided for a region may differ according to the GCM selected [32]. Some of them even provide contradictory information [33,34]. Therefore, ten climate projections have been selected. These projections come from five GCMs and two greenhouse gas (GHG) emissions scenarios, considered as representative concentration pathways (RCPs) [35] in the IPCC Fifth Assessment Report (AR5).
In this context, it is essential to quantify the alterations of the river caused by climate change both in magnitude and temporal variability. In fact, the choice of the appropriate hydrological indicators of river ecosystems may offer a robust tool in water management related to riparian habitats' long-term conservation. To date, more than 170 hydrologic metrics to describe aspects of flow regime and ecological status [36][37][38] have been published, although some redundancy is involved in the statistical formulation of some of them [37]. The indicators of hydrologic alteration (IHA) [8] are some of the most commonly used sets worldwide [39][40][41][42][43]. Indeed, Laize et al. [44] combined the use of the continental-scale water model WaterGAP (University of Kassel, Kassel, Germany) with an IHA-based methodology to quantify ecological risk due to the climate change flow alteration at the pan-European scale. Stagl and Hattermann [45] assessed the regions at the highest risk of climate change-driven flow alterations in the Danube river by using the Soil and Water Integrated Model [46] and eight ecological river flow indicators. Similarly to the IHA approach, Chatterjee et al. [47] evaluated climate change impacts on hydrologic flow regimes in the Great Plains of Kansas, combining the Soil and Water Assessment Tool (SWAT) and IHA. The European Union has also recommended [48] the use of the Indicators of Hydrologic Alteration in Rivers (IHARIS) [49,50]. This methodology has been applied both in Spanish and in Mediterranean river typologies inside and outside of Spain [51][52][53][54]. Both sets of metrics (IHA and IAHRIS) compare different flow regimes in spatial or temporal periods in order to assess variability in ecological parameters based on the inter-annual variation of hydrologic conditions due to the importance of flow in the health and diversity of aquatic ecosystems [7]. Both methodologies offer a software which, through the introduction of hydrologic series, obtains the results of the indicators. Though both metrics are adequate for the purpose of this study, the IAHRIS method has been selected [52,55] because it takes into account the unique Spanish climatic features, especially the pronounced inter-and intra-annual variability [49].
To our knowledge, although many studies related to the assessment of hydrological river alterations comparing land use changes [56][57][58], dam constructions [59][60][61], or climate change [44][45][46][47] have been developed worldwide, none of them have taken into account the range of alterations in ecological and hydro-geomorphological parameters when different climate projections are used, combining the Soil and Water Assessment Tool (SWAT) and IAHRIS. Thus, the goals of this research can be divided into three stages: (1) to compare the climatic variations in precipitation and temperature in short-, medium-and long-term period in northern Spain, (2) to assess the performance of the SWAT model in daily time interval and to generate future streamflow series according to the ten climate change projections, and (3) to analyze and discuss the range of hydrological alterations in the river selected using IAHRIS methodology. The contents of the paper are structured as follows: the study area is introduced in Section 2; the datasets used and methodology are described in Section 3; Section 4 presents the results and discussion; and Section 5 highlights the main conclusions.

Study Area
The Eo River is located in the north of Spain ( Figure 1). The watershed area is 819 km 2 and the length is 79 km, running into the Cantabrian Sea. The average annual streamflow is 18 m 3 /s, varying daily between 425 and 0.60 m 3 /s in 1971-2013 period. The Eo watershed lies mainly in a mountainous terrain whose altitudes are higher towards the north-south direction. The altitudes range between sea level and 1132 m above sea level (MASL), with an average of 473 m. Nevertheless, 48% of the area ranges from between 200 and 600 m in altitude and only the 0.5% is above 1000 MASL. Mean slopes are around 32%, so therefore plains are scarce and are mainly situated in the downstream river banks and some peaks are highly affected by erosion. The Eo watershed is characterized by the humid oceanic climate, with average temperatures around 13-14 °C and 1300 mm of mean yearly rainfall. The nearer to the mountainous areas, the higher the precipitation, where there is the influence of Cantabrian winds. Both human and physical factors of the watershed favor land use where hardwood forests and other forests account for more than 50% of the total watershed. Urbanized, industrial or anthropogenic surfaces do not reach 1% of the total area, and therefore traditional agricultural economy is usually preserved. Low population density and the conservation of traditional agrarian environments have turned 90% of the watershed into a Biosphere Natural Reserve, so land use change is not considered in this study because it seems unlikely that unsustainable development would be allowed in this protected area. Moreover, nearly 8% is qualified as a Community Interest Area (CIA) and the wetlands of the estuaries are resting places for many of the migrating birds, which are regarded as Special Protection Areas (SPAs) under the 'Birds Directive'. Likewise, the species of conservation concern are mainly constituted by taxa linked to river systems, indicating the key role of river ecosystems in the biodiversity preservations in the territory occupied by this Biosphere Natural Reserve. Concerning the modification of the river course and the quality of the body of water, the river basin organization considers the Eo River as 'natural' nature and 'good´ status, respectively [62].

SWAT Model
SWAT is a semi-distributed hydrological model [63]. SWAT divides a watershed into different sub-watersheds to take into account the river network and topography heterogeneity of the former. Likewise, the sub-watersheds are divided into hydrologic response units (HRUs), characterized by unique soil, land cover and slope class combination. In recent years, SWAT has been widely tested in numerous watersheds worldwide and has proven to have a satisfactory ability to model the water cycle [64,65] under very different climatic conditions and multiple land uses. However, the main weaknesses with SWAT are related to the nonspatial reference of the HRUs and the significant conceptual limitations in how groundwater flow and storage are simulated. Several authors have highlighted the simplistic representation of groundwater processes [66,67], which may cause poor model performance in groundwater-dominated watersheds. This is not the case in the Eo River basin, so the use of the SWAT model was considered appropriate.
The water balance equation (Equation (1)) used by SWAT considers as inputs daily precipitation and maximum and minimum temperatures.
where SWt is the final soil water content (mm), SWinit is the initial soil water content (mm), t is the time in days, Rday(i) is the precipitation on day i (mm), Qsurf(i) is the surface runoff (mm), Ea(i) is the evapotranspiration (mm), Wseep(i) is the percolation (mm), and Qgw(i) is the amount of baseflow (mm).

Model Set Up and Data Sets
For this study, 42 years of data (1971-2013) were used as the baseline period. Daily total precipitation (mm) were obtained from a gridded dataset from the Spanish State Meteorology Agency (AEMET), with a resolution of 5 km. Maximum and minimum temperatures (°C) were obtained from the gridded dataset called SPAIN02 and the AA-3D as interpolation method [68], with a resolution of 0.1°. The discharge data at the catchment outlet of this natural flow regime river were extracted from the Hydrographical Study Centre website [69].
The open source QGIS (QGIS Development Team, http://qgis.org) interface for SWAT, QSWAT, was used in this study [70]. SWAT requires hydro-meteorological, topography, soil properties, and land-use/land-cover inputs in the catchment. Heargaves' method [71] was used for simulating the potential evapotranspiration. Delineation of the Eo watershed and its stream network was based on a digital elevation model (DEM) obtained from the National Center for Geographic Information (IGN) with a 25 m grid cell resolution. Following the Food and Agriculture Organization of the United Nations (FAO) criteria for a watershed characterization [72], slopes were divided into three groups (<8%, 8-30%, >30%). Land cover was extracted from the European project Corine Land Cover 2012 and soil data from the Harmonized World Soil Database. The predominant soil types are Leptosols (60.5%). Other soil types present are Cambisols (37%) and Podzols (2.4%). The subdivision into HRUs was performed using threshold levels of 10% in soils, slopes and land uses, resulting in 7 sub-basins, including 288 HRUs.

Sensitivity Analysis, Calibration and Validation
Sensitivity analysis and calibration of parameters of the SWAT model were automatically carried out in SWAT-CUP (Texas A&M University, Texas, United States) using the SUFI-2 algorithm [73]. The selected parameters were calibrated using the observed daily discharge. Firstly, we did 500 model runs to obtain the most sensitive parameters. Afterwards, 2 iterations of 1000 simulations were run, readjusting the parameters after the second iteration. The input data series was divided into three phases: warm up, calibration, and validation. The model was calibrated daily for a period of seventeen years (1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990), using the Nash-Sutcliffe efficiency as an objective function [74]. After calibration, the model was validated using the daily discharge data of twenty-three years . Three years (1971)(1972)(1973) were used to warm up the model.

Goodness-of-Fit Tests
The results obtained both in calibration and validation phases were evaluated statistically using the widely used statistics Nash-Sutcliffe efficiency coefficient (NSE), percent bias (PBIAS), R 2 and root mean squared error (RMSE), defined in Table 1.

Performance Metric Equation Range
Optimal Oi is the ith observed data, � is the mean of the observed data, Ei is the ith estimated data, � is the mean of the estimated data and n is the total number of observations. Furthermore, following Kalin et al. [75], a grading method was established to evaluate the good performance of the model obtained based on the NSE and PBIAS values, as shown in Table 2.

Climate Change Scenarios
In order to simulate future scenarios, general circulation models (GCMs) are used to assess the possible changes to the climate in the future by modelling the effect of particular future emission scenarios of greenhouse gases on the climate [76]. As many studies have recommended using ensemble results of several climate change models [77], we have used the Climate Change Toolkit (CCT) software [78] (University of New South Wales, Kensington, Australia) which links with a global database of five GCMs from ISI-MIP (Inter-Sectoral Impact Model Intercomparison Project) [79], covering the 1950-2006 period for historical data and 2006-2099 for future data: GFDL-ESM2M (GCM1), HadGEM2-ES (GCM2), IPSL-CM5A-LR (GCM3), MIROC-ESM-CHEM (GCM4) and NorESM1-M (GCM5). Two representative concentration pathways (RCP4.5 and RCP8.5) have been selected for each GCM. Thus, ten climate change scenarios were used in the present study. All climate databases are at 0.5° spatial resolution.
Furthermore, a comparative analysis of the downscaled bias-corrected GCMs was carried out according to Pulido-Velazquez et al. [80] and Blanco-Gómez et al. [81], to compare the historical series from the climate stations with those from the GCMs. The indicator Id (Equation (2)) was used to assess the goodness-of-fit of every GCM, using the mean and standard deviation of the climate variables.
where subindex i indicates a specific GCM, V is the climate variable (V1 is precipitation and V2 is temperature), S1 is the monthly mean value, S2 the monthly standard deviation and j is the months for a mean year. [51]. It is composed of 25 parameters (Table 3); the first nine are obtained from the annual and monthly streamflow, respectively, while the rest are calculated with the daily flowrates. The results are classified as usual values, maximum extreme values (floods) and minimum extreme values (droughts). Furthermore, in each of these levels, four issues of high environmental significance are taken into account: 1) magnitude, which determines water availability, 2) frequency of an event, indicating the variability of a flowrate regime and a constraint of ecological and geo-morphological dynamics, 3) duration of an event, which is, in extreme situations, directly linked to the resilience of the species, 4) seasonality or regularity, related to the life cycles of the species. Using these parameters, a set of indicators (IAH) are formulated to compare the differences between natural and altered regimes ( Table 4) that, in this study, correspond to historic and predicted regimes, respectively. Mostly alteration indicators are calculated as the ratio between altered and natural regimes parameters values. They vary in the range of 0 to 1, where 1 means the unaltered state and 0 represents total alteration. Detailed description of the calculation method and the IAH expressions are presented in Appendix A. Five hydrologic statuses were established (Table 5) according to the range of indicator values and using the code of colours recommended for the ecological quality ratios (EQR) [82]. The lower the value of the index is, the higher the hydrologic alteration is and the status is very deficient. The opposite is also true: if the value of the index is higher, this indicates lower hydrologic alteration and excellent status. Table 5. Allocation categories criteria for partial alteration indicators (IAH).

Level IV Deficient
The IAHRIS software provides two types of result for each of the three components of the regime (usual, droughts and floods). On one side, there is a spider chart which shows both historic and predicted indices for the issue assessed, comparing the similarities or differences in indicators ( Figure  2). On the other side, a global index of alteration (IAG) is calculated for usual, drought and flood values. This global index is computed as the ratio between the areas of the historic and predicted polygons of the spider chart. Similarly, these indices of global alteration (IAG) are classified according a colour code (Table 6).

Meteorological Variations under Climate Change Scenarios
We analysed the changes in precipitation and both maximum and minimum temperatures of historical recorded data in the 1970-1999 period compared to the ten climate change scenarios studied (Table 7) according to three different time periods: short-term (2010-2039), mid-term (2040-2069) and long-term (2070-2099). Overall, the longer the period is, the less the average annual precipitation in nearly all the scenarios is, up to −22% in long-term for GCM1 RCP8.5. Notwithstanding, GCM3 shows barely any variation in RCP4.5 and, even in RCP8.5, an increase of 10% and 6% for the short-term and mid-term, respectively. Furthermore, trends of reduction in precipitation are also smoothed out in GCM2, especially in RCP4.5, where there is an increase of 10% in the short-term and the decrease in the long-term is only −1% compared to the historical period 1970-1999. Likewise, RCP 4.5 usually showed less variation compared to RCP8.5 in terms of total rainfall. Average annual rainfall in RCP4.5 is between 4 and 9% higher than RCP8.5 in the whole simulated period. GCM3 is the only model that presents different behaviour from the rest and average total precipitation in RCP8.5 increases 1% compared to RCP4.5, although the long-term value differs from previous periods, falling up to -14%. These uncertainties in the changes in precipitation were also shown by Ribalaygua et al. [83] in Aragón (northeastern Spain). Regarding average temperatures, all models but GCM1 predicted high increases as time goes on. These changes also depend on the scenario and whether minimum or maximum temperatures are considered. As with precipitation, both minimum and maximum temperatures rise higher as the following period is analysed. The average maximum temperatures are between +1.15 and +2.06 °C higher in the short-term compared to the historical period  and are very similar for RCP4.5 and RCP8.5, while they reach up to 1.81-4.21 °C and 3.21-7.01 °C in the long-term for RCP4.5 and RCP8.5, respectively. Table 7. Comparison of historical and predicted values under climate change of average annual precipitation (Pav), average maximum temperature (Tmax) and minimum temperature (Tmin). (ST, MT and LT refer to short-term, mid-term and long-term period, respectively).

GCM RCP Period Pav (mm) ∆P (%) Tmax (°C) ∆Tmax (°C) Tmin (°C) ∆Tmin (%)
Historic  1316. 30  Furthermore, a comparative analysis was carried out to evaluate the GCMs that fit best for historical data series (precipitation and temperature) from climate stations, according to the guidelines of Pulido-Velazquez et al. [80]. Table 8 shows the Id indicators for each GCM. The GCMs with the highest Id values had the worst fit. Considering the total Id value, GCM3 was identified as the best one, and GCM2 was the worst one. As the single Id values for precipitation and temperatures were similar (Table 8), higher differences were found in the standard deviations of temperature. Despite this, all of the GCMs were taken into account in this study because of the similarity of their total Ids. Finally, in order to analyse the range of annual average differences in models with respect to the historic period 1970-1999, Figure 3 plots the extreme scenarios (maximum and minimum deviations) for precipitation and maximum and minimum temperatures. Figure 3a shows that the least rainfall scenario (GCM1 RCP8.5 2070-2099 period) is defined by a nearly constant slope with an absolute minimum in the summer, as historical values. However, current peaks in spring and autumn seasons are smoothed. In contrast, precipitation in GCM2 RCP4.5 in the 2010-2039 period (short-term) is lower than the average historical values from June to November, but rainfall in the rest of the months is higher, especially in February and December. In spite of this, total precipitation in this scenario is 3% lower than historical data. Regarding temperatures, tendencies in both extreme scenarios are quite similar and the highest increases occur in the summer, as with the historical period. Nevertheless, the parallelism found between the hottest scenario and the historical one in maximum temperatures remains almost constant throughout the year with increases of +4.5°, except in the summer, when it may get as high as +10 °C in the long term. The model with the best fit (GCM3) also shows higher variations in both minimum and maximum temperatures, and all of the GCMs were also considered in this study as references for extreme climate conditions.

Sensitivity Analysis, Calibration and Validation of the SWAT Model
The first step in calibration was to identify the parameters that have greater influence in outflow results. Thus, a preliminary sensitivity analysis (500 model runs) was conducted with 20 parameters (Table 9), according to their occurrence on the objective function in the calibration step for the variable flow rate [84][85][86][87]. Afterwards, two iterations of 1000 simulations were run for the most sensitive parameters in the SWAT-CUP programme to make an auto-calibration using the SUFI2 algorithm. A ranking of parameter sensitivities was assessed with the p-value and the t-stat index ( Table 9). The former provides the significance of the sensitivity and the latter a measure of sensitivity. The most sensitive parameters are those that present the lowest value of p-value and the highest value of t-stat [72]. The ranking of the most sensitive parameters observed in this study (Table 9) was also supported by the previous findings in the north of Spain by Jimeno-Saéz et al. [84] Malagó et al. [88] and Raposo et al., [89]. Finally, the set parameters (Table 10) were also selected as the most relevant in other studies [89][90][91][92]. Evaluation of the parameters in Table 10 indicates that the SCS runoff curve number (CN2) decreased by 11% compared to the default value, thus increasing the infiltration and decreasing the runoff. Therefore, the values of ALPHA_BF and ALPHA_BNK obtained (0.87 and 0.91, respectively) indicate a slow response of runoff to rainfall due to the importance of aquifers in the study area. The calibrated values of groundwater parameters (GWQMN, GW_DELAY and REVAPMN) were similar to those used by Malagó et al. [88] for northern Spanish watersheds. The high value of Mannings' "n" value for overland flow (OV_N) is related to the abundant vegetation which is prevalent in northern Spain [89].

Model Performance
The values of NSE on a daily basis for calibration and validation were 0.73 and 0.62, respectively, which classifies the model as good according to the criteria listed in Table 2, since the PBIAS values were in both cases less than 25%, specifically 3.52% and 6.62% for the calibration and validation phases, respectively, indicating a slight but nearly negligible tendency to consider larger values than the observed ones. The other statistical indices also showed a good performance both in calibration and validation. R 2 was 0.74/0.65 and RMSE 10.73/14.71 for the calibration/validation stages. All these goodness-of-fit measures support the use of the model in predicting streamflow under the different climate scenarios previously described.

Variations in Streamflow
Following an analysis and comparison of the results of the total streamflow in the different scenarios studied (Table 11), all of them show a decrease with respect to the historical period, and the differences range between −5% and −53%. The lowest variation is achieved with GCM2-RCP4.5-ST because it is one of the scenarios with a higher increase in precipitation. However, the highest decrease is also shown for the same GCM2, but in RCP8.5 and long-term, since it is one of the scenarios with the combination of both less precipitation and a higher increase in average temperature. Although this scenario had the worst fit among the selected scenarios, it showed similar behaviour compared to GCM3 (best fit) because the variations ranged from −12% to −39% in RCP8.5 and from −22% to −25% in RCP4.5. In general, as climatic values showed, the later the period is, the lower the level that the streamflow reached. Similar conclusions are usually shown between RCP 4.5 and 8.5. Although the disparities between scenarios are evident, the average decrease in the GCMs predicts a reduction of around 20-30% of water resources in the current century, as has also been suggested both for Europe [93][94][95][96] and for Spain [97]. Table 11. Comparison between historical and predicted average annual streamflow (Q). (ST, MT and LT refers to short-term, mid-term and long-term period, respectively). According to the annual average decrease of around 20-30% of streamflow in the 21st century, the flow-duration curve of GCM3-RCP4.5-MT (∆Q = −22%) is analysed in Figure 4 in order to assess which types of streamflow are the most affected in climate change scenarios. The major difference between the historic regime (1970-1999) and GCM3-RCP4.5-MT appears for percentiles lower than 25%, which indicates that peak flows will be severely reduced and, consequently, so will parameters related to flooding, such as periodic flood plain events and thus, the development and the variety of species, especially at their earliest stages. On the contrary, the percentiles higher than 25% are fairly similar, and therefore the impact on both usual and drought values is expected to be weak.  Table 12 shows IAG results and the spider-charts for usual values and extreme cases, both in floods and droughts, according to the IAHRIS method. In general, IAG values range between 0.75 and 0.37 for usual values, 0.33 and 0.08 for floods and 0.50 and 0.28 for droughts. This indicates that the extreme values in floods are more deficient than both the usual and extreme values in droughts. In fact, nearly 25% of the future climate scenarios assessed showed an excellent status for usual values, obtaining an average of 0.57. In the case of droughts, the average is reduced to 0.38, which means a good status, though the threshold for moderate is 0.36 and therefore 37% of the scenarios are classified as moderate in droughts. On the other hand, the average of the floods values is 0.17, close to deficient status. In fact, the deficient status is obtained in floods in 60% of cases and none of the other scenarios reach the good status. On the contrary, for both usual values and droughts, IAG is never below 0.28, resulting in a moderate status in the most pessimistic scenario. These results of usual values are indicative of the general availability of water and the pattern of temporal variability, which means that the quality and heterogeneity of habitat will be able to regulate the intrusion and expansion of exotic species [98]. Furthermore, moderate alteration in droughts ensures habitat availability and the sustainability of hydro-habitat in dry periods due to the maintenance of water levels and velocities and quality of water [99][100][101][102].

IAHRIS
When taking into account the considered period, indices values often worsen in droughts and usual values the later the period is, due to the decrease in the streamflow, although this variation between short-term and long-term in studied scenarios is, on average, below the 20% threshold both for the usual and drought values. By contrast, the index improves slightly for the floods or remains invariable the later the period is.
The spider-charts support the interpretation of previous findings depending on the displacement of the vertices of future scenario from historic regime figures, used as reference. For the usual values, the most unfavourable parameters are the magnitude of monthly volume and annual variability, which are found to be moderate and, although annual and monthly streamflow are not directly concerned with any particular geomorphological or ecological stage, monthly volumes and inter-annual variability of streamflow are clearly decisive to the ecosystem [103,104]. Furthermore, the reduction in most of the monthly usual streamflow values implies a decrease in the available hydro-habitat [105], which will be sharpened at the end of the century. With respect to the rest of the parameters of usual values, greater closeness between polygons is shown in the short-term, with values of IAH between 0.47 and 0.75. However, the inner polygon shrinks the later the period is, especially with regard to monthly magnitudes, but scarcely in annual ones. Indeed, the inter-annual disparities of these values do not often exceed 40%, whereas, when taking into account monthly differences, they increase up to 60%, indicating that, while it is true that the annual volumes have been moderately reduced, they considerably fall in the monthly period, especially in quantity compared to the historic period of reference, as is shown in Figure 4. On the contrary, the seasonality of maximums and minimums (alteration indicators E1 and E2) reaches the value of 1 in most cases, which means excellent conditions regarding this issue. This synchronicity sets the pace of vital processes and are linked to a set of environmental factors that are essential for the temporal maintenance of diversity of habitats as well as the stimulation of germination and dispersal [106]. In general, predictions in usual values confirm minor to slight decreases in all scenarios both in magnitude and variability. Moreover, extreme variability, characterized as the difference between the maximum and minimum input that has occurred in the year, does not show serious alterations, as is also reflected in Figure 4.
Concerning flood situations, all GCMs show deficient status, except in GCM2 and GCM3, where the mean value of 0.26 corresponds to moderate status. The poor results achieved for all scenarios are mainly due to the values of zero both in the frequency of connectivity flow (IAH9) for all scenarios and the duration of floods (IAH13) for GCM1, GCM4 and GCM5, which is in line with previous studies [28][29][30]. This connectivity ensures aquatic access to the riparian and flood plain river system and the maintenance of suitable moisture conditions for different biota growth stages [107], also stimulating riparian habitat rejuvenation, succession processes in the riparian forest [108,109], meander sequences, and lateral channel flushing [110], etc. Furthermore, the loss of the connectivity with flood plains implies a continuous aging of riparian habitat, endangering the renewal of the species [111]. Nevertheless, these findings must be interpreted with caution given the uncertainty of SWAT models in estimating the large values of streamflow [84,112]. Concerning other magnitudes of flood indicators (IAH7, IAH8 and IAH10), a decrease with respect to the historic regime is also observed, which is more pronounced in IAH7 and IAH8 than in IAH10. This indicates that all geomorphological functions related to the effective discharge (QGL) will be affected and to the cleaning and the regeneration substrate functions of usual floods, though to a lesser extent [113]. This change in magnitude will be accentuated with the loss of diversity of these events and is reflected in the IAH11 and IAH12 indicators in all scenarios. The seasonality of floods (IAH14) remains unaltered compared to the historic regime, while their durations (IAH13) show different results depending on the scenario considered: highly altered in GCMs 1, 4 and 5 and slightly changed in GCMs 2 and 3.
Regarding droughts, the number of days without zero flow in the river (IAH20) always reaches a value of 1, reflecting the continuous flow in the Eo river for all the projected scenarios. Furthermore, flow related to droughts (QS and Q95%) decreases with regard to the historical regime, as is shown in the IAH15 and IAH16 indicators. It is also observed that the duration of droughts (IAH19) is severely increased for all the scenarios, reaching extreme values of consecutive dry days the longer the period taken into account is. This could affect the resilience capacity of species [114]. Likewise, GCM1, GCM4 and GCM5 show moderate modifications in variability of usual droughts for all periods. These facts endanger the availability of habitat, as well as the characteristics of hydro-habitat, such as loss of habitat for aquatic organisms, loss of stream connectivity, as well as deterioration of water quality, alteration of food resources, and changes in interspecific interactions [115,116].
Despite the excellent-good status in usual values and good-moderate status in drought values, the poor results in flood situations would not contribute to the cleaning and revitalization of the substrate, which could inhibit the diversity of the hydro-habitat, as well as reduce the stimulus for the migratory movements that ensure the accessibility to breeding areas [117][118][119]. It should also be mentioned that, although much of the flow quantification can be assessed in linear terms, the ecosystems' dynamics are mainly ecologically driven [120]. Table 12. IAG and spider-chart for the combination GCM-RCP-Period (GCM-RCP-P).

Conclusions
Climate change is an indisputable reality all over the world. This change has the potential to result in an important decrease in river flows, both their magnitude and variability. The deviation from the natural regime will affect the river ecosystem and the availability and diversity of aquatic habitats. This study aims to assess the hydrologic alteration caused by different climatic change scenarios. Ten climate scenarios were taken into account from five GCM and two RCP for each one. The results of meteorological parameters in most scenarios in the studied area pointed out a reduction in rainfall (up to −22%) and an increase in temperatures, both maximum (1-7 °C) and minimum ones (1-4 °C). Overall, the longer the period is, the more important are the differences compared to the historic regime.
The SWAT model was applied on a daily basis to a northern Spain watershed in order to provide a hydrological model to be used in obtaining the projected future flow rates. The statistical results of calibration are NSE = 0.73, PBIAS = 3.52, R 2 = 0.74 and RMSE = 10.73. The validation results are NSE = 0.62, PBIAS = 6.62, R 2 = 0.65 and RMSE = 14.71. These results are indicative of the SWAT model's good performance. After the calibration phase, predicted streamflow values were obtained for the ten climate change scenarios, resulting in an incremental decrease the longer the term considered, varying between −5% (in short-term) and −53% (in long-term) depending on the scenario assessed.
IAHRIS shows ( Figure 5) that alteration for usual and drought values remains between excellent (23%) and good (77%) and good (63%) and moderate (37%) status, respectively, becoming usually worse the longer the term is and linked to predicted meteorological parameters. Furthermore, the increase in the duration of the droughts may endanger availability of habitat and water quality during these periods. Notwithstanding, these relatively soft climatic changes in usual and dry conditions are determinant in the maintenance of quality and diversity of habitats and in the survival of the most sensitive species in critical environmental periods. Furthermore, moderate alterations help to govern the intrusion of alien species. However, flood values show a moderately deficient regime (40-60%, respectively), which implies an instability of river morphology, a progressive reduction in the section of the river and an advance of the aging of riparian habitat, endangering the renewal of the species. Thus, in spite of the results for usual and drought values under the studied climate scenarios, this Biosphere Natural Reserve runs the risk of the loss of the diversity of the hydrohabitat and the loss of the connectivity with the flood plains due to the reduction in the magnitude of floods and, consequently, the cleaning and revitalization of the substrate where the river ecosystem is developed. This study intends to focus on the consequences of climatic change in river ecosystems through the joint use of SWAT and IAHRIS. The assessment and quantification of the alterations in usual and extreme values reveal that the maintenance of the natural regime involves not only usual values, but also extreme streamflow that may affect the quality and diversity of riparian habitats and geomorphology of the natural course of the rivers. This methodology can be applied in any watershed to warn of the potential environmental risks of climatic change and to adopt specific measures that could prevent the endangerment of the hydro-habitat. The effects of groundwater, land use change, drainage and harvesting rates, for example, could also be assessed. However, because they do not have an important role, they were not included in the modelling, although they can be considered in future research or for more accurate analysis.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A. IAHRIS Methodology
The indicators of alteration (IAH) are calculated as the quotient between the value of the characterization parameter in an altered regime and the same parameter's value in a natural regime, as shown in Equation (A1).
The IAH indicators used in the current study are defined and developed below. M1 (Magnitude of the annual volumes): Evaluates the distortion caused by the circulation to the magnitude of the annual discharges, in comparison with their value in a natural regime, obtained by Equation (A2): where ����� is the average annual volume in an altered regime and ����� is the average annual volume in a natural regime.
M2 (Magnitude of the monthly volumes): In order to avoid potential compensations between different months of the year when working with annual values, this indicator is proposed to mitigate the former state by allocating the same weight to the alteration in each month. The alteration is therefore assessed month by month, for each month of each year. The indicator is finally expressed as the mean value of all of the estimated monthly alterations (Equation (A3)): where ( ℎ ) ������������������� is the average monthly volume of month i in an altered regime and ( ℎ ) ������������������� is the average monthly volume of month i in a natural regime.
M3month (Magnitude of the monthly volumes, from each month): The alteration is assessed individually for each month to provide the 12 monthly indicators, as indicated in Equation (A3): V1 (Variability of the annual volumes): Evaluates alterations in variability in terms of annual discharges by using Equation (A4): where CV(AA)n is the coefficient of variation of the series of annual volumes in an altered regime and CV(AA)n is the coefficient of variation of the series of annual volumes in a natural regime. V2 (Variability of the monthly volumes): Evaluates changes in variability in terms of monthly discharges (Equation (A5)): where ( ℎ ) is the coefficient of variation of the series of monthly volumes for month i in an altered regime and ( ℎ ) is the coefficient among the variation of the series of monthly volumes for month i in a natural regime.
V3 monthi (Variability of the monthly volumes, for each month): Alteration is assessed individually for each month to provide the 12 monthly indicators by using Equation (A6): V4 (Extreme variability): Evaluates distortions caused by the circulating regime to the extreme variability of habitual values in comparison with a natural regime, calculated using Equation (A7): where Am maximum, year i, a is the maximum monthly volume of year i in altered regime (and similarly, the minimum for Am minimum year i, a.), Am maximum, year i, n is the maximum monthly volume of year p in a natural regime (similarly for Am minimum, year p, n, N is the number of years available in an altered regime and M is the number of years available in a natural regime. E1 (Seasonality of maximum values): Evaluates the distortion caused by the circulating regime to the seasonality of the monthly maximum natural discharges. After the table is built of relative frequencies for the presence of maxima each month, the change in the maxima's seasonality is estimated as the delay or difference (measured in months) between the month with the highest maximum frequency in a natural regime and the month with the maximum frequency in an altered one, as given by Equation (A8): E2 (Seasonality of minimum values): Evaluates the distortion caused by the circulating regime to the seasonality of the monthly natural minimums. After the table is built of relative frequencies for the presence of minima each month, the change in the minima's seasonality is estimated as the delay or difference (measured in months) between the month with the highest maximum frequency in a natural regime and the month with the maximum frequency in an altered one, as given by Equation (A9): IAH7 (Magnitude of the maximum floods). Evaluates changes in the average value of circulating floods. The index assesses changes in the magnitude of peak floods, in which the parameter is the average of the maximum annual flow series in each type of regime, by using Equation (A10): where ����� is the mean of the maximum annual flows within the available series in an altered regime and ����� is the mean of the maximum annual flows within the available series in a natural regime. IAH8 (Magnitude of the effective discharge): Evaluates the change to the values of flows with special geomorphological significance. The indicator defined in Equation (A11) uses the effective discharge as its reference parameter and corrects the relationship yielded by the exponent 0.5, taking into account the direct relationship between the stream width and the square root of QGL: where QGLa is the effective discharge corresponding to an altered regime and QGLn is the effective discharge corresponding to a natural regime. IAH9 (Frequency of the connectivity discharge): Evaluates the change in the frequency of flows that ensure periodic connection with the floodplain, as indicated in Equation (A11): where QCONEC is the connectivity discharge in a natural regime, equivalent to the flow that duplicates the effective discharge's return period in this regime, Ta(QCONEC) is the return period corresponding to connectivity discharge in an altered regime and Tn(QCONEC) is the return period corresponding to connectivity discharge in a natural regime. The Gumbel extreme value distribution is considered to characterize the largest flows.  where Q5% a is the 5% exceedance percentile on the flow-duration curve in an altered regime and Q5% n is the 5% exceedance percentile on the flow duration curve in a natural regime.
IAH11 (Variability of maximum floods): Quantifies the distortion in the interannual variability of maximum flows. This index uses the coefficient of variation among the maximum annual flows (Qc) as its estimator, as given by Equation (A12): where ( ) is the coefficient of variation for the series of maximum annual flows in an altered regime and ( ) is the coefficient of variation for the series of maximum annual flows in a natural regime.
IAH12 (Variability of the flushing floods): Quantifies the distortions in the interannual variability of flushing floods (Q5%). This indicator (Equation (A13)) uses the coefficient of variation of the Q5% series of values in a natural and an altered regime as its estimator.
where ( 5% ) is the coefficient of variation of the series of values corresponding to a flushing flood in an altered regime (annual series of Q5%) and ( 5% ) is the coefficient of variation of the series of values corresponding to a flushing flood in a natural regime (annual series of Q5%).
IAH13 (Flood duration): Ascertains whether flood periods maintain similar durations in altered regimes to durations in natural regimes. This index is defined as the maximum number of consecutive days when the (Q5%)n threshold is equalled or surpassed. Changes to flood duration are consequently evaluated as the quotient between the value of this parameter in an altered regime and the corresponding figure in a natural regime, as shown in Equation (A14): where (Maximum consecutive days when Q > Q5%n)a is the maximum number of consecutive days, on average, in an altered regime, when a flushing natural flood is equalled or surpassed and (Maximum consecutive days when Q > Q5%n)n is the maximum number of consecutive days, on average, in a natural regime, when a flushing natural flood is equalled or surpassed. IAH14 (Flood seasonality): This indicator is proposed for evaluating changes in the seasonal aspect of floods. The procedure focuses on assessing the effects on flood seasonality for each month of year, and then calculating the change indicator as the average of these monthly indicators. The monthly indicator IAH 14 month i is obtained as follows: if the absolute value of the difference between a natural regime and an altered regime is higher than 5, then IAH 14 month i = 0, otherwise, it is calculated as given by Equation (A15): where ABS is the absolute value, NATURAL is the average days per month when Q > Q5% n in a natural regime for the month considered and ALTERED is the average days per month when Q > Q5% n in an altered regime for the month considered. Indicator IAH 14 is provided, which is generated as the average of the 12 monthly values. IAH15 (Magnitude of the extreme drought): Evaluates changes in the minimum annual flows. The indicator estimates changes in drought magnitude on the basis of the relationship between the mean minimum annual flow corresponding to altered and natural one, as given by Equation (A16).
where ����� is the mean minimum annual flow of the available series in an altered regime and ����� is the mean minimum annual flow of the available series in a natural regime. IAH16 (Magnitude of the habitual drought): Evaluates changes in drought magnitudes by using Equation (A17), while avoiding the sole use of extreme minimum values, given that they camouflage other situations that, while possibly less critical, are nevertheless of great environmental importance: where 95% is the flow (m 3 /s) corresponding to a habitual drought in an altered regime and 95% is the flow (m 3 /s) corresponding to a habitual drought in a natural regime. In both cases, the normal drought is defined by the flow value corresponding to the 95% percentile in the respective duration-flow curve, which is similar, in terms of duration, to the value for flows that, on average, are not equaled or surpassed on 18 days per year.
IAH17 (Variability of the extreme drought): Evaluates changes in the variability of minimum annual flows. The parameter given by Equation (A18) is used to estimate the coefficient of variation of the minimum annual flows corresponding to each of the regimes: where ( ) is the coefficient of variation of the minimum annual flows corresponding to an altered regime and ( ) is the coefficient of variation of the minimum annual flows corresponding to a natural regime.
IAH18 (Variability of the habitual drought): Evaluates changes in the variability of the most frequent droughts. The indicator (Equation (A19)) uses the coefficient of variation of the series of values Q95% corresponding to the two considered regimes: where ( 95% ) is the coefficient of variation of the annual series of values corresponding to a habitual drought in an altered regime (annual series of Q95%) and ( 95% ) is the coefficient of variation of the annual series of values corresponding to a habitual drought in a natural regime (annual series of Q95%).
IAH19 (Drought duration): Tests whether drought periods in an altered regime have a similar duration to those in a natural one by using Equation (A20): where ( ℎ < 95% ) is the maximum number of consecutive days in an altered regime, on average, when a normal natural drought is not surpassed, and ( ℎ < 95% ) is the maximum number of consecutive days in a natural regime, on average, when a normal natural drought is not surpassed. IAH20 (Number of days with null flow): Evaluates changes in the characteristic number of days with null flow in a natural regime. The procedure centres on evaluating effects on the duration of noflow periods for each month of the year and then calculating the alteration indicator as an average of these monthly indicators. The monthly indicator IAH20 month i is obtained as follows. If the absolute value of the difference between a natural and an altered regime is higher than 5, then IAH20 month i = 0; otherwise, it is calculated as given by Equation (A21): where ABS is the absolute value, NATURAL is the average days per month when Q = 0 in a natural regime for the month under consideration and ALTERED is the average days per month when Q = 0 in an altered regime for the month under consideration. IAH21 (Drought seasonality): Evaluates changes in drought seasonality. As with IAH20, the reference flow chosen to determine the presence or absence of a drought in a given month is Q95% in a natural regime. The procedure centres on evaluating the effects on the seasonality of droughts for each month of the year, and then calculating the change indicator as an average of these monthly indicators. The monthly indicator IAH21 month i is obtained as follows. If the absolute value of the difference between a natural and an altered is higher than 5, then IAH21 month i = 0, otherwise, it is calculated as given by Equation (A22).
where ABS is the absolute value, NATURAL is the average days per month when Q < Q95% n in a natural regime for the considered month and ALTERED is the average days per month when Q < Q95%