Waning habitats due to climate change: effects of streamflow and temperature changes at the rear edge of the distribution of a cold- water fish

Climate change affects aquatic ecosystems altering temperature and precipitation patterns, and the rear edge of the distribution of cold-water species is especially sensitive to them. The main goal was to predict in detail how change in air temperature and precipitation will affect streamflow, the thermal habitat of a cold-water fish (brown trout, Salmo trutta Linnaeus 1758), and their synergistic relationships at the rear edge of its natural distribution. 31 sites in 14 mountain rivers 15 and streams were studied in Central Spain. Models at several sites were built using regression trees for streamflow, and a non-linear regression method for stream temperature. Nine global climate models simulations for the RCP4.5 and RCP8.5 (Representative Concentration Pathways) scenarios were downscaled to a local level. Significant streamflow reductions were predicted in all basins (max. -49 %) by the year 2099, showing seasonal differences between them. The stream temperature models showed relationships between models parameters, geology and hydrologic responses. Temperature was sensitive to 20 the streamflow in one set of streams, and summer reductions contributed to additional stream temperature increases (max. 3.6oC), although the most deep-aquifer dependent sites better resisted warming. The predicted increase in water temperature reached up to 4.0oC. Temperature and streamflow changes will cause a shift of the rear edge of the species distribution. However, geology conditioned the extent of this shift. Approaches like these should be useful in planning the prevention and mitigation of negative effects of climate change by differentiating areas based on the risk level and viability of fish 25 populations.

in fish energy balance, affecting the rate of food intake, metabolic rate and growth performance (Forseth et al., 2009;Elliott and Elliott, 2010;Elliott and Allonby, 2013). It is also involved in many other physiological functions such as blood and reproductive maturation (Jeffries et al., 2012), reproductive timing (Warren et al., 2012), gametogenesis (Lahnsteiner and Leitner, 2013), cardiac function (Vornanen et al., 2014), gene expression (White et al., 2012;Meshcheryakova et al., 2016), ecological relationships (Hein et al., 2013;Fey and Herren, 2014), and fish behaviour (Colchen et al., 2016). 5 Natural patterns of water temperature and streamflow are profoundly linked with climatic variables (Caissie, 2006;Webb et al., 2008). Stream temperature is strongly correlated to air temperature (Mohseni and Stefan, 1999), but streamflow has a complex relationship with precipitation (McCuen, 1998;Gordon et al., 2004). In addition, temperature influences the precipitation type (rainfall or snowfall) and snowmelt occurrence, and inversely, river discharge is also a main explanatory factor of water temperature for some river systems (Neumann et al., 2003;van Vliet et al., 2011;Toffolon and Piccolroaz, 10 2015). Furthermore, geology affects surface water temperature by means of the groundwater discharges (Caissie, 2006), influencing the aquifer deep (shallow or deep) and the water's residence time (Kurylyk et al., 2013, Snyder et al., 2015. Climate change is already affecting aquatic ecosystems altering temperature and precipitation patterns. Stream temperature increases have been documented for the last decades throughout the globe, in Europe (e.g., Orr et al., 2015), Asia (e.g., Chen et al., 2016), America (e.g., Kaushal et al., 2010) and Australia (e.g., Chessman, 2009). Flow regimes are also being 15 influenced by changes in precipitation, although those trends vary by climatic region (IPCC, 2013;Morán-Tejeda et al., 2014). The predictions of the International Panel of Climate Change (IPCC, 2013) suggest that these alterations will continue throughout the XXI century and it will have consequences on the freshwater fish distribution (e.g., Comte et al., 2013;Ruiz-Navarro et al., 2016). This alteration may especially affect cold-water fish, which have been shown to be very sensitive to climate warming (Williams et al., 2015;Santiago et al., 2016). For example, among salmonids, DeWeber and Wagner, 20 (2015) found temperature to be the most important determinant of the occurrence probability of the brook trout, Salvelinus fontinalis (Mitchill, 1814).
The rear edge of the distribution (sensu Hampe and Petit, 2005) of a cold-water species is especially sensitive to changes in temperature in addition to reductions in the available habitable volume (i.e., streamflow). The impact of water temperature on the distribution of salmonid fishes is well documented (e.g., Beer and Anderson, 2013;Eby et al., 2014); however, the 25 combined effects of temperature rise and streamflow decay remain relatively unexamined; with some exceptions, e.g. Wenger et al., (2011) or Muñoz-Mas et al., (2016). Jonsson and Jonsson (2009) predicted that expected effects of climate change on water temperatures and streamflow will have implications for migration, ontogeny, growth and life-history traits of Atlantic salmon, Salmo salar Linnaeus, 1758, and brown trout, Salmo trutta Linnaeus, 1758. Thus, investigation of these habitat variables in the context of several climate scenarios should help scientists to assess the magnitude of these changes 30 on the suitable range and life history.
The objective of this study was to predict how and to what extent the availability of suitable habitat for the brown trout, a sensitive cold-water species, will change within its current natural distribution under the new climate scenarios through the study of the changes in streamflow and temperature and their interaction. To this end hydrologic simulation with M5 models trees coupled with non-linear water temperature models at the daily time step were feed with a fine grained downscaling of the air temperature and precipitation predicted for the most recent climate change scenarios (IPCC, 2013). It is hypothesized that the flow regime and geological nature of the basin influences the thermal regime of rivers and streams; thus, the effects of basin geology on the temperature models and on the changes in thermal regimes were studied. Finally, the alterations of trout thermal habitat were assessed by studying the violation of the tolerance temperature thresholds of the brown trout. 5

Materials and methods
The logical framework followed is summarized in Fig. 1. First, more recent daily general climate models at local resolutions were downscaled. Then, the obtained local climate models were applied to create streamflow and temperature models. The results are daily values that might be useful for the assessment of hydrological resources and the changes in water quality and aquatic communities based on the changes in resources availability. As noted, this study is focused on fish habitat 10 suitability and availability. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License.

Study sites
A total of 31 sites at 14 mountain rivers and streams inhabited by brown trout was chosen with the aim of encompassing a diverse array of geological and hydrological conditions in the centre of Spain (between 39º53' N and 41º21' N latitudes): the Tormes River and its tributaries the Barbellido River, Gredos Gorge and the Aravalle River (Duero basin), Cega River, Pirón River (Pirón is a tributary of the Cega River, in the larger Duero basin), Lozoya River, Tagus River, Gallo River, Cabrillas 5 (Tagus basin), Ebrón River, Vallanca River (Vallanca is a tributary of the Ebrón River, in the Turia basin), Palancia River and Villahermosa River (Fig. 2). Major components of geology lithologically characterize the mountain sites in the Duero and Lozoya basins as igneous, altitudinally lower sites in the Duero basin are detrital (Cenozoic), and the eastern basins (Tagus, Gallo, Cabrillas, Ebrón, Vallanca, Palancia and Villahermosa) are characterized as mainly carbonated (Mesozoic).

15
Land use is mainly forest in all study site basins (CORINE Land Cover 2006[European Environmental Agency, 2007).
Only the lower basins of the downstream sites on the Cega and Pirón rivers are mosaics of forest and croplands, whereas the uppermost sites within the Tormes River basin (Barbellido and Gredos Gorge) lie above the current tree line. Land use change scenarios were assumed to not be probable. The studied reaches are not effectively regulated (only small weirs or natural obstacles exist). One large dam lies on the Pirón River (Torrecaballeros Dam; capacity: 0.324 hm 3 ; maximum depth: 20 26 m; and altitude: 1390 m a.s.l.), but it does not significantly alter the temporal streamflow pattern (Santiago et al., 2013).
Hydrological data characterized the streamflow regimes as extreme winter/early spring (groups 13 and 14 in the classification of Haines et al., [1988]). However, hydrographs show a west to east smoothing gradient (Fig. 3). Smoothing is associated with carbonated geology, whereas greater seasonality is associated with igneous and detrital geology.

Data collection 5
In each study site, water temperature was registered every two hours using 31 Hobo® Water Temperature Pro v2 (Onset®) and Vemco® Minilog data loggers located at several sites along the studied rivers (Table 1). Meteorological data were obtained from 9 thermometric and 15 pluviometric stations of the Spanish Meteorological Agency (AEMET) network, and 10 gauging stations (from the Water Administrations' official network) and were considered to model running flows. The closest AEMET-thermometric stations to the stream temperature monitoring sites with at least 30 years between 1955 and 10 the present were selected. The selected pluviometric stations were those located within the target river basin and upstream of the corresponding gauging station (Table 2).

Climate change modelling and downscaling
Data from nine global climate models associated with the 5 th Coupled Model Intercomparison Project were used, namely: BCC-CSM1-1; CanESM2; CNRM-CM5; GFDL -ESM2 M; HADGEM2-CC; MIROC-ESM-CHEM; MPI-ESM-MR; MRI-15 CGCM3; NorESM1-M (Santiago et al., 2016) which provided daily data to simulate future climate change corresponding to the Representative Concentration Pathways RCP4.5 (stable scenario) and RCP8.5 (more increasing CO 2 scenario) established in Taylor et al. (2009) and used in the 5 th Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2013). An array of nine general climate models was used to avoid bias due to the particular assumptions and features of each model, (Kurylyk et al., 2013). Historical simulations (XX century) were used to control the quality of the procedure 20 and to compare the magnitude of the predicted changes. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. Pourmokhtarian et al. (2016) remark the importance of the use of fine downscaling techniques; thus, a two-step analogue statistical method (Ribalaygua et al., 2013) was used to downscale the daily climatic data: the maximum and minimum air temperatures and the precipitation for each station and for each day. For both air temperature and precipitation, the procedure begins with an analogue stratification (Zorita and von Storch, 1999) in which the n days most similar to each problem day to be downscaled are selected using four different meteorological large-scale fields as predictors: (1) the speed and (2) direction 5 of the geostrophic wind at 1000 hPa and (3) the speed and (4) direction of the geostrophic wind at 500 hPa. In a second step, the temperature determination was obtained using a multiple linear regression analysis using the selected n of most analogous days. This was performed for each station and for each problem day, as well as for maximum and minimum temperatures. The linear regression uses forward and backward stepwise selections of the predictors to select only the predictive variables for that particular case. For precipitation, a group of m problem days (the whole days of a month were 10 used) were downscaled together: 'preliminary precipitation amount' averaging the rain amount of its n most analogous days was obtained for each problem day. Thus, the m problem days from the highest to the lowest 'preliminary precipitation amount' could be sorted. For assigning the final amount of rain, each one of the amounts of rain of the m × n analogous days was taken. Then those m × n amounts of rain were sorted, and then those amounts were clustered in m groups. Every quantity was then assigned, orderly, to the m days previously sorted by the 'preliminary precipitation amount'. Further 15 details of the methodology are described in Ribalaygua et al. (2013).
A systematic error is obtained when comparing the simulated data from the climate models with the observed data, inherently associated with every downscaling methodology and every climate model (both of which usually introduce bias into the data). To correct this systematic error, the future climate projections were corrected according to a parametric quantile-quantile method (Monjo et al., 2014), performed by comparing the observed and simulated Empirical Cumulative 20 Distribution Functions (ECDF) and linking them by the ECDF of the downscaled European Centre for Medium-Range Weather Forecasts ERA-40 reanalysis daily data (Uppala et al., 2005).
As a result, for each climate change scenario, daily maximum and minimum temperatures (used to infer the mean temperature) and precipitation were obtained for each climate model which were used as input to simulate the runoff and water temperature under these climate change scenarios. 25

Hydrological modelling
The prospective prediction of the runoff flows was performed with data-driven hydrological models developed with the M5 algorithm (Quinlan, 1992). M5 has been demonstrated to be proficient to model daily streamflow (Solomatine and Dulal, 2003;Taghi Sattari et al., 2013), even in studies involving climate change (Muñoz-Mas et al., 2016). The M5 algorithm is a tree-like recursive partitioning technique, similar to Classification and Regression Trees (CARTs) (Breiman et al., 1984) or 30 Random Forests (RF) (Breiman, 2001). Recursive partition iteratively divides the input space into smaller subspaces, usually maximising the homogeneity of the child nodes. This partitioning typically begins with the most discriminant split and proceeds to the least one and is applied in a hierarchical fashion to each of the new branches of the tree until the maximum Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. number of allowed partitions or any other constraint is achieved (Vezza et al. 2015). However, instead of assigning a single value or category to each terminal nodes (i.e., leaf), M5 fits a multi-linear regression model (Quinlan 1992). Consequently, the ultimate tree becomes a piecewise multi-linear model, which can be seen as a committee of linear models with each member specialized in particular subsets of the input space (Taghi Sattari et al., 2013). Based on the multi-linear models at the leaves, M5 allows extrapolation, in contrast with other machine learning techniques that have demonstrated a poor or null 5 capacity to do so (Hettiarachchi et al., 2005).
The M5-based hydrological models were developed in R (R Core Team, 2015) with the Cubist package (Kuhn et al., 2014). One single M model tree was trained for each gauging station and predictions were supported by the nearest observation (i.e., neighbours = 1) to reduce the potential prediction of unreliable flows. Finally, M5 was allowed to decide the ultimate number of leaves in the model tree (i.e., the number of models in which the input space is eventually divided). 10 Following previous studies, the hydrological models were trained using daily, monthly and quarterly data lags of historical rainfall and temperature data as input variables (Solomatine and Dulal, 2003;Taghi Sattari et al., 2013;Muñoz-Mas et al., 2016). These three groups of variables were intended to reflect the causes of peak, normal and base flows. The study encompassed several rivers that may present different hydrologic behaviours; therefore, the starting set of input variables was larger than that used in other studies (Solomatine and Dulal, 2003;Taghi Sattari et al., 2013;Muñoz-Mas et al., 2016). 15 Daily variables included the precipitation and temperature from the current day to the 15 th previous day (16 variables in total). The monthly variables were calculated using the moving average for the 12 previous months (12 variables in total), and the quarterly data were calculated from the moving average for the current month to the 24 th previous month (8 variables in total). Consequently, daily variables overlapped with the current month variable, and the first four quarterly variables overlapped with the monthly data. Thus, 72 variables were used (36 each for temperature and precipitation). 20 The whole set of input variables may be relevant for some river systems, although it may cause M5 to overfit the data in others. Therefore, unlike previous studies that used smaller input variable sets (Muñoz-Mas et al., 2016), the ultimate variable subset was optimized following the forward stepwise approach (Kittler, 1978). This greedy approach relies on iteratively adding inputs while performance on the test data set improves and stopping as soon the performance degrades.
However, this approach may cause consideration of unrelated variable sets (i.e., disjoint rainfall and temperature variable 25 periods). To address such potential inconsistency, the optimization began by testing the rainfall-related variables and only tested the temperature variables for lags coinciding with those rainfall-related variables selected. No precaution was taken regarding correlation among inputs (Solomatine personal communication), and the forward stepwise approach sought the maximization of the Nash-Sutcliffe Efficiency (NSE) index (feasible range from −∞ to 1; Nash and Sutcliffe [1970]) in a fivefold cross-validation as indicated in Borra and Di Ciaccio (2010). Following previous studies (Fukuda et al., 2013;30 Muñoz-Mas et al., 2016), once the optimal variable set was determined, single M5 models trees were developed to perform the prospective prediction of the flow rate under the climate change scenarios.
To assess the significance of the streamflow trends throughout the century, Sen's slope was used (Trend package of R (Pohlert, 2016); p-value ≤ 0.05) with horizons H-2050 and H-2099.
Finally, the variation of the patterns of the monthly mean streamflow was studied by means of cluster analysis of the change rate of the normalized monthly mean streamflow at H-2050 and H-2099, for the RCP4.5 and RCP8.5 scenarios (Ward Hierarchical Clustering -cluster package of R [Maechler, 2013]). 5

Stream temperature modelling
Daily Mean Stream Temperature (DMST) was used because it better reflects the average conditions that fish (particularly trout) will experience for an extended period of time (Santiago et al., 2016) and also it averaged daily fluctuation in the radiation and heat fluxes.
A non-linear regressive model was used to determine DMST from Daily Mean Air Temperature (DMAT). The Mohseni et 10 al. (1998) regressive logistic model was preferred to machine-learning models such as classification and regression trees (De'ath and Fabricius, 2000) and random forests (Breiman, 2001) because parameters in the former have physic significance and higher transferability (Wenger and Olden, 2012). Thus, models described temperature behaviour adjusting a sigmoid curve to the data affecting the trend of the higher water temperatures, which present an attenuated gradient due to the evaporative water-cooling, being more reliable than linear regressions (Mohseni and Stefan, 1999). At the same time, 15 prediction on air temperature will overpass the actual temperature range thus models which assure transferability as nonmachine learning should be preferred (Wenger and Olden, 2012). The Mohseni et al. (1998) model was modified by Santiago et al. (2016) to include the trajectory of temperature to improve its performance when daily temperature is modelled (Term 1 in Eq. [1]). In the modified model Ts is DMST, Ta is DMAT, µ is the minimum stream temperature (given in ºC), α is the maximum stream temperature (ºC), β represents air temperature at which the rate of change of the 20 stream temperature with respect to air temperature is maximum (ºC), γ (ºC -1 ) is the value of the rate of change at β, and λ is a coefficient (dimensionless) representing the resistance of DMST to change with respect to 1-day variation in DMAT (ΔTa).
Negative values of λ are due to the resistance to temperature change and, thus, must be subtracted from the expected temperature: the more resistant the stream is to temperature change, the closer λ is to zero. The less resistant to change, the more negative λ is. 25 In addition, the model was further modified by introducing the daily mean flow (Q) through a logistic function with three parameters (Term 2 in Eq. [1]): ω is the maximum observable variation in temperature due to the flow difference (given in 30 ºC), τ represents the flow value at which the rate of change of the stream temperature with respect to the flow is maximum (m 3 ·s -1 ), and δ (m -3 ·s) is this maximum rate at τ.
A blockwise non-parametric bootstrap regression (Liu and Singh, 1992) was used to estimate the parameters of both the modified Mohseni models (with and without streamflow), and residual normality and non-autocorrelation were checked with Shapiro test and Durbin-Watson test, as well the 7 day-lag PACF (Partial Autocorrelation Function) were obtained.
Calculations were performed using R. A 95 % confidence interval was calculated for each parameter. Bayesian information criterion (BIC) and Akaike information criterion (AIC) were used to test the eight-parameter models against the five-5 parameter models.
The parameter µ was allowed to be less than zero in the modelling process even though this is the freezing temperature.
Thus, the function would truncate at the freezing point. The relationship between the thermal amplitude α-µ and the indicator of thermal stability λ was studied using the Pearson correlation.
Geology is determinant of the residence time of deep-groundwater in the aquifer (Chilton, 1996), and residence time 10 influences the discharge temperature. To explore the relationships between thermal regimes and geology, a stratified study of both the geology classes of the parameter values was completed by means of a t-test with Bonferroni correction (p-value < 0.05). In the same sense, increments of the annual averages of the daily mean (ΔTmean), minimum (ΔT min ) and maximum (ΔT max ) stream temperatures were calculated and studied by geological classes.
The variation of the patterns of the monthly mean stream temperature was studied by means of cluster analysis of the 15 temperature increments at H-2050 and H-2099, for the RCP4.5 and RCP8.5 scenarios (Ward Hierarchical Clusteringcluster package of R [Maechler, 2013]).

Thermal habitat changes
Several tolerance temperatures have been described for brown trout (Elliott et al., 1995;Forseth et al., 2009;Elliott and Elliott, 2010). The elected threshold for this study was the occurrence of DMST above 18.7ºC for seven or more consecutive 20 days because it has proven to be the most realistic value to represent the realized thermal niche (Santiago et al., 2016). The minimum period of 7 consecutive days is usually the established time for determining thermal tolerance (Elliott and Elliott, 2010), and when this period is exceeded, the death risk increases harshly. The use daily mean rather than daily maximum is better to reflect the average conditions that trout must experience for an extended period. In addition, studying DMST for 7 consecutive days above the threshold because is the more accurate option (Santiago et al., 2016). 25 Once DMST was modelled, the frequency of events of 7 or more consecutive days above the threshold per year (Times To assess the general trend in thermal habitat alterations at the middle (H-2050) and the end of the century (H-2099), the TAT≥7, DAT and MCDAT were calculated at each sampling site for each climate change scenario and compared with 30 current conditions. In this way, the sampling sites spread to a wide range of geological and hydrological conditions in three main basins. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License.
The sampling sites distribution at Cega, Pirón and Lozoya streams allowed us to study the relation between the annual average of DMST and altitude, and a high correlation was detected (R 2 = 0.986, 0.985 and 0.881, respectively). A digital elevation model (at a 5-m resolution, obtained from LIDAR, IGN [National Geographic Institute of the Spanish Government]) was used to perform an altitudinal interpolation of the model parameters to determine the water temperature along the stream continuum to simulate the climate change scenarios' effects, and then to obtain the percentage of 5 stream/river length lost for trout. ArcGis ® 10.1 (by ESRI ® ) was used to manage the DEM.

Climate change
Under the climate change scenarios, all the meteorological stations experience important increases in temperature (DMAT) through the century, and this trend proved to be more pronounced in the RCP8.5 scenario. The main changes will occur in 10 summer, and to a lesser extent, in winter (annual trends in Fig. 4; seasonal results per location in Fig. S1 to S24 [Supplementary Material 1]). The temperature variations run parallel in both scenarios until mid-century, when the RCP8.5 scenario predict a similar trend and the increases become smaller in the RCP4.5 scenario: the annual change in temperatures for the RCP4.5 fluctuates between 2ºC and 2.5ºC at mid-century (3-4ºC in summer) and between 2.5ºC and 3.5ºC at the end of the century (3.5-4.5ºC in summer), and the annual change for the RCP8.5 is between 2ºC and 3ºC at mid-century (3. [5][6][7][8][9][10][11][12][13][14][15] 4.5ºC in summer) and between 5ºC and 7ºC at the end of the century (7-8ºC in summer).
In absolute terms (mm·day -1 ), annual precipitation fluctuates around the no-change value (zero) (Fig. 4), but seasonal trends differ (Fig. S1). The changes in precipitation are mainly small, but large changes are detected at some sites, especially in autumn, with important precipitation decreases in the RCP8.5 scenario. For the same time scale, RCP4.5 returns to the current values at the end of the century. In relative terms (%), the changes are also small but sufficiently significant to be 20 considered: there is an annual change in the RCP4.5 from 0 to -15 % at mid-century (0 to -20 % for the October-March period) and 10 to -10 % at the end of the century (10 to -10 % for the October-March period), and there is an annual change in the RCP8.5 from 10 to -10 % at mid-century (0 to -20 % for the October-March period) and 0 to -20 % at the end of the century (0 to -20 % for the October-March period).
Daily mean air temperatures of the ensemble models for each meteorological station are shown in the Supplementary 25 Material 2 (Dataset S1).

Hydrological regimes
The hydrological models performed well-all of them achieved NSE values above 0.7 with the exception of Tormes-Hoyos (NSE=0.673); nonetheless, each model was constituted by different variable sets (Table S1 [

Geographical pattern
The current streamflow patterns, strongly modulated by geology, will experience a progressive increase of the influence of 10 climate change that is reflected in the differences among the clustered sites (Fig. 6).
In the RCP4.5 scenario at H-2099 the clustering is very weak (agglomerative coefficient: 0.56) and heterogeneous. At H-2050, the clustering yields a geographical change pattern; however, in spite of the distance, Ebrón station is aggregated to Tormes stations because the similarity of the hydrological response (a.c.: 0.73).
The agglomerative coefficient (0.61) for clustering in the RCP8.5 scenario at H-2050 is weak, with the only recognizable 15 pattern grouping the stations that will lose more flow throughout the year. The most important changes will be produced in the RCP8.5 scenario at H-2099 (a.c.: 0.72). Cega-Lastras and Lozoya form the group of the most heavily affected sites by flow reductions; one of them is in the Duero basin, and the other in the Tagus basin, but they are geographically not far from one other. Another group share with the former the generalized losses excepting at summer. The remaining group can be disaggregated in two other: one with the stations located at Tormes basin and Ebrón which show significant losses at spring 20 and autumn, and increases at summer, and the second group showing homogeneous but no so heavy losses. In this cluster the geologic component is attenuated.
In general, streamflow losses will increase throughout the century but not consistently in all the sites. At the most western (Tormes) and eastern (Ebrón) stations, flow values will increase again at H-2099 after decreasing in the fifties. Lozoya will suffer the most substantial flow losses, followed by Pirón and Cega-Lastras, Tagus and Gallo, and Cabrillas. 25

Model parameter behaviour and general trends
In accordance with the BIC and the AIC calculated for each temperature model, the inclusion of the running flow improves model performance in 12 out of 31 study cases. In the other 19 cases, the five parameters model was used (Table 4). In cases in which no improvement occurs with the model that included streamflow effect, at times, no convergence of values is 30 observed in the regression process, and at other times, the obtained values do not improve the results as the streamflow Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License.   Table 5, and daily mean stream temperatures from the climate change models are in the Supplementary Material 5 (Dataset S3).
As the thermal amplitude increased, λ become more negative (less resistance). As the thermal amplitude decreases, λ 5 approaches zero (more resistance).
By the end of the XXI century, the predicted average increase in the mean annual stream temperature among the sites is 1.1ºC for the RCP4.5 scenario simulations (range 0.3-1.6ºC) and 2.7ºC for the RCP8.5 (range 0.8-4.0ºC). Annual maximum average temperature increases are predicted to be 0.8ºC for the RCP4.5 (range 0.1-1.5ºC) and 1.6ºC for the RCP8.5 (range 0.2-3.0ºC), and annual minimum average temperature increases are 1.0ºC (range 0.4-1.8ºC) and 2.7ºC (range 1.1-4.5ºC), 10 respectively. The most important increases are predicted to occur in winter, with summer experiencing lesser increases.
In RCP4.5 scenario, ΔT min at quaternary-detrital sites behaves significantly different than in carbonated and igneous sites. 5 This difference is only found between quaternary-detrital and carbonated sites in RCP8.5 scenario. ΔT mean exhibits significant differences between the carbonated and igneous and quaternary detrital sites, in both scenarios. All of these results are common at H-2050 and H-2099. Finally, at H-2050, significant differences are found between ΔT max in carbonated and igneous sites in both scenarios and, while this occurs in the same way at H-2099 in RCP4.5 scenario, the differences are also significant between carbonated and quaternary-detrital sites in RCP8.5 scenario (Fig. 8). 10

15
The results of the cluster analysis of the monthly mean stream temperature of sites show highly homogeneous groups between the different combinations of horizons and scenarios, being the thermal response of the rivers and streams very linked to geology. The carbonated sites from Cabrillas stream to the East and Pirón 3 (under the strong influence of a calcareous spring) form a group of sites that shows low thermal amplitude and in which lambda is closest to zero. In the other extreme, a group mainly formed by granitic sites (Lozoya and Tormes basins, and various sites in the detrital basin of 20 Cega-Pirón basin) shows higher thermal amplitude and lower values of lambda than the former group. All other sites have intermediate values of both, thermal amplitude and resistance.

Effect of streamflow reduction on stream temperature
The effect of streamflow change on temperature is analysed at the following sites: Tormes 2, Tormes 3, Pirón 1, Cega 1, Lozoya 1 to 4, Cabrillas, Ebrón 1 and Vallanca 1 and 2. These are the sites in which the 8-parameter model improves upon 5 the 5-parameter model. In all cases, stream temperature differences between the five and eight parameter models are found, and summer flow reduction led to an increase in stream temperature, increasing DAT, TAT≥7 and MCDAT. Among these cases, the threshold is only surpassed in Lozoya and Tormes, increasing the thermal habitat loss. In Cega 1, Cabrillas and Ebrón, α is below the thermal threshold, and in Pirón 1, the temperature increase is not sufficient to exceed the threshold.
The differences between the annual maximum stream temperatures for both models at H-2050 and H-2099 imply increases 10 up to 3.6ºC in Tormes 2, with smaller increases in Tormes 3 (2.5ºC) and the Lozoya sites (0.9ºC to 1.7ºC). Because of the improved model, increases in the maxima of 1.1ºC in Cega and 1.4ºC in Pirón are observed. At the Ebrón, Vallanca and Cabrillas sites, the increases in the maxima are smaller (0.01ºC to 0.3ºC).

Effect of climate change on brown trout thermal habitat
In the predictions at H-2050, the 18.7ºC threshold (TAT≥7) is violated at 8 sites in the RCP4.5 scenario and 6 sites in the 15 RCP8.5 scenario. In the climate change simulations for the end of the XXI century (2099 horizon), the threshold is violated at 8 sites in the RCP4.5 scenario and 13 sites in the RCP8.5 scenario. These violations imply the loss of thermal habitat for brown trout.
By the end of the century (H-2099), the most significant increases in TAT≥7 (Fig. 9)  habitat, affecting up to 56 %, 11 % and 66 % of the stream length, respectively. In the Cega and Pirón rivers, the loss is referred the total main channel length inhabited by the trout (98 and 77 km at Cega and Pirón, respectively). In Lozoya, the Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. length lost is the reach (20 km) located upstream of the Pinilla reservoir, which disconnects the reaches of the Lozoya upstream and downstream of the dam. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License.

Climate change
Results from the 5 th Assessment Report of the Intergovernmental Panel on Climate Change and its annexed Atlas of Global and Regional Climate Projections (IPCC, 2013) for the Mediterranean area suggest that droughts are not likely to increase in 5 the near future. However, air temperature is expected to rise, subsequently increasing evapotranspiration. Our downscaled results predict greater temperature increments than did the IPCC (2013), which can produce more relevant ecological consequences (Magnuson and Destasio, 1997;Angilletta, 2009). In the Mediterranean area, regional studies have used coarser resolution, which may be appropriate for their goals (e.g., Thuiller et al., 2006). However, they can fail when more local predictions are needed. Thus, fine downscaling techniques as applied in this study must be used when a finer prediction 10 resolution is needed.

Streamflow
Worldwide, abundant information is available regarding the impact of recent climate change on streamflow regimes (e.g., Luce and Holden 2009;Leppi et al., 2012) and, more specifically, in our study area (e.g., Ceballos-Barbancho et al., 2007;Lorenzo-Lacruz et al., 2012;Morán-Tejeda et al., 2014), but detailed predictions are uncommon (e.g., Thodsen, 2007). At 15 the regional level, a reduction in water resources is expected at the Mediterranean area (IPCC, 2013). Milly et al. (2005) predicted a 10-30 % decrease in runoff in Southern Europe in 2050. In another global scale study, van Vliet et al. (2013) predicted a decrease in the mean flows of greater than 25 % in the Iberian Peninsula area by the end the century (2071-2100), using averages for both the SRES A2 and B1 scenarios (Nakicenovic et al., 2000). Our results predict that mean flow will be similar to that value (-23 %, range: 0-49 %), although the emissions scenarios in this study are more severe (more increasing 20 CO2 scenario) than in previous studies.
The predictions for the RCP4.5 scenario show changes ranging from non-significant to significant flow reductions of up to 17 %. In the RCP8.5 scenario, significant predictions of reduction were more generalized, ranging up to 49 % of the annual streamflow losses. Our results also suggest a relevant increment in the number of days with zero-flow for some stations in the detrital area in this scenario (RCP8.5). The most sensitive variable was Q min , which showed significant declines at a 25 larger number of sites. The predicted changes are compatible with those obtained in previous studies, although they were done for a greater scale (as cited: Milly et al., 2005;van Vliet et al., 2013). Apparent differences between the streamflow reductions in this study and the former (with lower flow reductions in the present study) might be because the former predictions were for the whole region (Iberian Peninsula), whereas ours are focused in mountain reaches; the scale of focus is likely also important. First, the climate and then geologic features are determinant for the hydrologic response (Raghunath, 30 Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Brown trout is a sensitive species to changes in discharge patterns because high intensity floods during the incubation and emergence periods may limit recruitment (Lobón-Cerviá and Rincón, 2004;Junker et al., 2014). In the Iberian Peninsula, extreme discharges during winter are expected to increase (Rojas et al., 2012) and, therefore, negatively affect trout 5 recruitment. Thus, the predicted changes in the hydrological regime can subject brown trout populations to greater stochasticity, which may occasionally drive some populations to insuperable bottlenecks. Trout is a polytypic species and might have adaptable phenology and rather high intra-population variability in its life history traits that might allow a resilient response to variation in habitat features (Gortázar et al., 2007;Larios-López et al., 2015), especially in the marginal ranges (Ayllón et al., 2016). However, despite these strong evolutionary responses, the current combination of warming and 10 streamflow reduction scenarios is likely to exceed the capacity of many populations to adapt to new conditions (Ayllón et al., 2016). Consistently with regional predictions (Rojas et al 2012;Garner et al., 2015), significant flow reductions are expected during summertime in most of the studied rivers and streams at the end of the century, and this may mean, in turn, the reduction of available habitat for trout (Muñoz-Mas et al., 2016). The increase of extreme droughts involving absolute water depletion in certain reaches of the streams may be critical for some trout populations. 15 Bustillo et al. (2013) recommended the assessment of the impacts of climate change on river temperatures using regressionbased methods that rely on logistic approximations of equilibrium temperatures (Edinger et al., 1968), which are at least as robust as the most refined classical heat balance models. Because our models are based on this argument, the stream temperature models developed in this study were found to be robust. 20

Stream temperature
However, we also looked for relationships between thermal regime and other environmental variables in addition to air temperature and streamflow, such as geology. Bogan et al. (2003) showed that water temperatures were uniquely controlled by climate in only 26 % of 596 studied stream reaches. Groundwater, wastewater and reservoir releases influenced water temperature in the remaining 74 % of the cases. Loinaz et al. (2013) quantified the influence of groundwater discharge in the temperature variations in the Silver Creek Basin (Idaho, USA), and they conclude that a 10 % reduction in the groundwater 25 flow can cause an increase of over 0.3ºC and 1.5ºC in the average and maximum temperatures, respectively. Our studied reaches did not suffer influences from wastewater or reservoir releases (with the exception of the Torrecaballeros Dam on the Pirón River). Kurylyk et al. (2015) showed that the temperature of shallow groundwater influences the thermal regimes of groundwater-dominated streams and rivers. Since groundwater is strongly influenced by geology, we can expect it to be a good indicator of the thermal response, as shown here. The parameters of the models accurately described the thermal 30 performance of the study sites, and we found significant relations among parameters, geology and hydrologic response.
Thermal amplitude (α-μ) and temperature at maximum change rate (β) were lower, and the resistance parameter (λ) was closer to zero, in basins with highly buffered aquifers (mainly carbonated) than in the others, especially in igneous basins. Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. Arismendi et al. (2014) concluded that regression models based on air temperature can be inadequate for projecting future stream temperatures because they are only surrogates for air temperature. In this study, we show that parameters of both the Mohseni model and our modification integrate the information of other variables, such as geology or flow regime. The parameter values could change with changes in geology (for instance), but these changes could be modelled by studying the behaviour of the parameters. The dynamics of parameters is a new research front, especially their asymptotic behaviour, 5 because these models may less accurately predict extreme temperature events (Arismendi et al., 2014). However, this may affect our study only if the value of the asymptote lies below the tested thresholds of fish tolerance. This can occur in highly buffered thermal regimes as the deep-groundwater very dependent streams. Thus, the limitation of ignoring the effect of climate warming on groundwater (subsurface water and deep water) must not be forgotten. Thermal sensitivity of shallow groundwater is different between the short-term (e.g., seasonal) and the long-term (e.g., multi-decadal), and the relationship 10 between air and water temperature is not necessarily representative of that. This variability must be taken into account to not underestimate the climate warming effects (Kurylyk et al., 2015). New research must be done to solve this limitation of the used model by integrate this criterion in the parametrization. Since a fine solution could need prohibitive methods (Kurylyk et al., 2015) losing the advantages that make attractive the model (input data easy to get), a compromise between improved precision and increased cost must be met. 15 Regressive models are substantially site specific compared to deterministic approaches (Arismendi et al., 2014). However, parameters of these regressive approaches are still physically meaningful, and need a lower amount of variables that can limit its feasibility in areas where data are scarce. Consequently, the value of this type of model is its applicability to a large number of sites where the only available data are for air temperature (and precipitation and streamflow to a lesser extent).
Our results show that predictions are improved when streamflow is included in the water temperature model, especially in 20 the periods when low flow and high temperatures are coincident, as is often true in temperate rivers and streams (Arismendi et al., 2012). Some rivers show no sensitivity to the introduction of streamflow in the model, which must not necessarily be due to the lack of influence of flow on the water temperature but rather to its minor relevance compared to other sources of noise.
In the group of streams that showed sensitivity to the introduction of streamflow into the temperature model, differences 25 between the stream temperatures estimated by the two models (with and without streamflow) were important in granitic basins, and summer flow reductions further increased stream temperatures. Other studies highlight the influence of streamflow on water temperature (e.g., regulated rivers: Sinokrot and Gulliver, 2000; groundwater discharge influence: Loinaz et al., 2013;and natural, regulated and snowmelt-fed rivers: Toffolon and Piccolroaz, 2015). When this occurs, peak flows must be modelled separately (Toth et al., 2000;Solomatine and Xue, 2004). Therefore, further development of our 30 model should include the effect of these rainfall peak flow events.
The predicted increase in water temperature will be important at most of the study sites. Annual mean variation rates will increase with time. Stewart et al. (2015) predict an increase by 1-2ºC by mid-century in the 80 % of stream lengths in Wisconsin and by 1-3ºC by late-century in 99 %, with a significant loss of suitable areas for cold-water fishes. The results of Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. Stewart et al. (2015) are not different from ours except that we expect greater increases at the end of the period (up to 4ºC).
Our results are also compatible with those by van Vliet et al. (2013), who predicted a temperature increase >2ºC in the Iberian Peninsula area by the end of the century; however, our results are more specific and precise. In this sense, Muñoz-Mas et al. (2016) also obtained similar results to ours in one river in Central Spain at the mid-century (daily mean flow reductions between 20-29 % and daily mean stream temperature increases up to 0.8ºC). However, we predict that minima are 5 more sensitive to climate warming than maxima.
The predicted increase in winter temperatures can affect the sessile phases of trout: i.e., eggs and larvae. These phases are very sensitive to temperature changes because it affects their physiology and because their durations are temperature dependent (e.g., Lobón-Cerviá and Mortensen, 2005;Lahnsteiner and Leitner, 2013). Thus, changes in the duration of the incubation and yolk sac absorption can affect the emergence time and, in turn, the sensitivity to hydrological regime 10 alterations (Sánchez-Hernández and Nunn, 2016). An increase in temperature can also reduce hatchling survival (Elliott and Elliott, 2010). In accordance with the herein presented results, the predicted synergy of streamflow reductions and water temperature increases will cause substantial losses of suitable fish habitat, especially for cold-water fishes such as brown trout (Muñoz-Mas et al., 2016).
The increases in threshold violations were important in our simulations. The duration of warm events (temperature above the 15 threshold) increased by up to three months at the end of the century in the most pessimistic scenario (RCP8.5). A continuous analysis of the whole-river response should be conducted to allow a spatially explicit prediction and to identify of reaches where thermal refugia are likely to remain. However, our results suggest that trout will not survive in these reaches because the existence of thermal refugia will be improbable or insufficient. In the Cega, Pirón and Lozoya rivers, losses of thermal habitat will be important and could jeopardize the viability of the trout population. Behavioural thermoregulatory tactics are 20 common in fish (Reynolds and Casterlin, 1979;Goyer et al., 2014); for instance, some species perform short excursions (<60 min in experiments with brook char, S. fontinalis) that could be a common thermoregulatory behaviour adopted by cold freshwater fish species to sustain their body temperature below a critical temperature threshold, enabling them to exploit resources in an unfavourable thermal environment (Pépino et al., 2015). Brown trout can use pool bottoms during daylight hours to avoid the warmer and less oxygenated surface waters in thermal refugia (Elliott, 2000). Nevertheless, if the warm 25 events became too long, the thermal refugia could become completely insufficient, thus compromising fish survival (Brewitt and Danner, 2014;Daigle, et al., 2014).

The brown trout distribution
According to our results, streamflow reductions synergistically contribute to the loss of thermal habitat by increasing daily mean stream temperature. It is especially relevant in summer in the Mediterranean area, when warmest temperatures and 30 minimum flows usually occur. The existence of thermal refugia is an alternative to fish survival and the probability for a water body to become a thermal refugium is highly geologically dependent. In our simulations, the most deep-aquifer Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License. dependent sites (Mesozoic-calcareous basins) will better resist warming. The habitat retraction at the rear edge of the actual distribution of brown trout is deduced to be geologically mediated.
The mountains of Central and Southeast Spain hold a rear edge of the brown trout native distribution (Kottelat and Freyhof, 2007). Fragmentation and disconnection of populations by newly formed thermal barriers may aggravate the already significant losses of thermal habitat by reducing the viability of populations and increasing the extinction risk. Thus, the rear 5 edge of the trout population in the Iberian Peninsula might be shifted to the northern mountains to varying extents depending on relevant mesological features, such as geology. The calcareous mountains of Northern Spain could be a refuge for trout by combining both issues: favourable geology and a relatively more humid climate. Caused by this differential response, the western portion of the Iberian range (plutonic, less buffered) will eventually experience more frequent local temperaturedriven extinction events, thus producing a greater shift northward, than in the Eastern Iberian end (calcareous, highly 10 buffered) of this range, which will remain more resilient to these local extinction events. However, the predicted streamflow reductions may act synergistically, reducing the physical space, and this may jeopardize the less thermally exposed populations. In the Iberian Peninsula, stream temperatures will raise less in the central and northern mountains than in the central plateau, and less in karstic than in granitic mountains. At the same time, the Mediterranean facade is expected to be more sensitive to warming and streamflow reductions than the Atlantic facade. Thus, brown trout populations in karstic 15 mountains in Northern Spain (Cantabrian Mountains and calcareous parts of Pyrenees) are better able to resist the climate warming than eastern populations in the granitic portion of Pyrenees (Santiago 2017). The same may occur in other parts of Southern Europe. Probably, the less pronounced thermal response of rivers and streams in the karstic areas will allow a better persistence of the brown trout population, although it would be also linked to changes in streamflow regimes.
Studying the great European basins, Lassalle and Rochard (2009) predicted brown trout "to lose all its suitable basins in the 20 southern part of its distribution area (Black Sea, the Mediterranean, the Iberian Peninsula and the South of France), but likely to continue being abundant in northern basins". Almodóvar et al. (2011) estimated that brown trout will be eradicated of almost the entire stream length of the studied basins in North Spain, and Filipe et al. (2013) estimated in 57 % the expected loss of the studied reaches at the Ebro basin (Northeast Spain). Our study shows important thermal habitat reduction, yet not so dramatic, of Iberian brown trout populations in mountain areas: the number of used General Climate Models, the 25 reliability of the downscaling procedure, the resolution of the stream temperature and streamflow models, and the method to study the threshold imply a substantial improvement of the detail (Santiago et al., 2016) over previous work. It is reasonable to infer that many mountain streams appear poised to be refugia for cold-water biodiversity this century (Isaak et al., 2016).
Approaches such as ours should be useful in planning the prevention and mitigation of the negative effects of climate change on freshwater fish species in the rear edge of their distribution. A differentiation of areas based on risk level and viability is 30 necessary to set management goals. Our results show that trout persistence requires a study of both temperature and streamflow dynamics at fine spatial and temporal scales. Managers need easy-to-use tools to simulate the expected impacts and the management options to address them, and the methods and results we provide could be a reliable way of doing so.  Table 1. Description of the data logger (thermograph) sites, specifying given name, UTM-coordinates (Europe WGS89), altitude (m above the sea level), code of the nearest temperature meteorological station with suitable time series for this study (AEMET: Spanish Meteorological Agency), orthogonal distance between the data logger and the meteorological station, number of recorded days for stream temperature and characteristic geological nature of the data logger site (the latter was obtained from ITGE, 2000). Bold letters indicate sites associated to the gauging stations.  Runoff annual -13.6 -15. 5 -25.8 -12.2 -12.7 Runoff annual -7.7 -12. 5 -12.5 -37.7 -48.6 -39 -32.1 -19 -17.5  Hydrol. Earth Syst. Sci. Discuss., doi: 10.5194/hess-2016-606, 2017 Manuscript under review for journal Hydrol. Earth Syst. Sci. Published: 9 January 2017 c Author(s) 2017. CC-BY 3.0 License.