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

Abstract. Climate changes affect aquatic ecosystems by altering temperatures and precipitation patterns, and the rear edges of the distributions of cold-water species are especially sensitive to these effects. The main goal of this study was to predict in detail how changes in air temperature and precipitation will affect streamflow, the thermal habitat of a cold-water fish (the brown trout, Salmo trutta), and the synergistic relationships among these variables at the rear edge of the natural distribution of brown trout. Thirty-one sites in 14 mountain rivers and streams were studied in central Spain. Models of streamflow were built for several of these sites using M5 model trees, and a non-linear regression method was used to estimate stream temperatures. Nine global climate models simulations for Representative Concentration Pathways RCP4.5 and RCP8.5 scenarios were downscaled to the local level. Significant reductions in streamflow were predicted to occur in all of the basins (max. −49 %) by the year 2099, and seasonal differences were noted between the basins. The stream temperature models showed relationships between the model parameters, geology and hydrologic responses. Temperature was sensitive to streamflow in one set of streams, and summer reductions in streamflow contributed to additional stream temperature increases (max. 3.6 °C), although the sites that are most dependent on deep aquifers will likely resist warming to a greater degree. The predicted increases in water temperatures were as high as 4.0 °C. Temperature and streamflow changes will cause a shift in the rear edge of the distribution of this species. However, geology will affect the extent of this shift. Approaches like the one used herein have proven to be useful in planning the prevention and mitigation of the negative effects of climate change by differentiating areas based on the risk level and viability of fish populations.


Introduction
Water temperatures are a primary influence on the physical, chemical and biological processes in rivers and streams (Caissie, 2006;Webb et al., 2008) and, subsequently, the organisms that live completely or partially in the water.Temperature is a major feature of the ecological niche of poikilothermic species (e.g.Magnuson and Destasio, 1997;Angilletta, 2009) and a key factor in energy balance of fish.It affects 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 function 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), Published by Copernicus Publications on behalf of the European Geosciences Union.ecological relationships (Hein et al., 2013;Fey and Herren, 2014), and fish behaviour (Colchen et al., 2017).
Natural patterns of water temperature and streamflow are profoundly linked with climatic variables (Caissie, 2006;Webb et al., 2008).Therefore, stream temperature is strongly correlated with air temperature (Mohseni and Stefan, 1999), whereas streamflow has a complex relationship with precipitation (McCuen, 1998;Gordon et al., 2004).In addition, atmospheric temperature influences the type of precipitation (rain or snow) that occurs and the occurrence of snowmelt; conversely, river discharge is also a main explanatory factor of water temperature for some river systems (Neumann et al., 2003;van Vliet et al., 2011).Furthermore, geology affects surface water temperatures by means of groundwater discharge (Caissie, 2006;Loinaz et al., 2013), influenced by the aquifer depth (shallow or deep) and the water's residence time (Kurylyk et al., 2013;Snyder et al., 2015).
Climate change is already affecting aquatic ecosystems by altering water temperatures and precipitation patterns.Stream temperature increases have been documented over the last several decades over the whole globe, such as in Europe (e.g.Orr et al., 2015, documented a mean increase in stream temperature of 0.03 • C per year in England and Wales), Asia (e.g.Chen et al., 2016, documented a mean increase in stream temperature of 0.029-0.046• C per year in the Yongan River, eastern China), America (e.g.Kaushal et al., 2010, documented mean increases in stream temperature of 0.009-0.077• C per year) and Australia (e.g.Chessman, 2009, documented mean increases in stream temperature of 0.12 • C per year between macroinvertebrate sampling campaigns).Abundant information is also available regarding the impact of recent climate changes on streamflow regimes worldwide (e.g.Luce and Holden, 2009;Leppi et al., 2012) and, more specifically, in the Iberian Peninsula (e.g.Ceballos-Barbancho et al., 2008;Lorenzo-Lacruz et al., 2012;Morán-Tejeda et al., 2014).However, detailed predictions are uncommon (e.g.Thodsen, 2007).The predictions of the Intergovernmental Panel on Climate Change (IPCC, 2013) suggest that these alterations will continue throughout the 21st century, and they will have consequences for the distribution of freshwater fish (e.g.Comte et al., 2013;Ruiz-Navarro et al., 2016).These changes may have an especially strong effect on 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 (2015) found stream temperature to be the most important determinant of the probability of occurrence of brook trout, Salvelinus fontinalis (Mitchill, 1814).
The rear edge populations (sensu Hampe and Petit, 2005: "populations residing at the current low-latitude margins of species' distribution ranges") of a cold-water species are especially sensitive to changes in water temperature, in addition to reductions in the available habitable volume (i.e.streamflow).The rear edge is the eroding margin of the range where lineages mix, the genetic drift and local adaptations increase, and droughts put populations under stress.The impact of water temperatures on the distribution of salmonid fish is well documented (e.g.Beer and Anderson, 2013;Eby et al., 2014); however, the combined effects of rising stream temperatures and reductions in streamflow remain relatively unexamined, with some exceptions (e.g.Wenger et al., 2011;Muñoz-Mas et al., 2016).Jonsson and Jonsson (2009) predicted that the expected effects of climate change on water temperatures and streamflow will have implications for the 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 on the suitable range and life history of these species.
The objective of this study is 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 a study of changes in streamflow and temperature and their interactions.Specifically, in this paper, we (i) assessed the effects of both streamflow and geology on stream temperature; (ii) predicted the changes in streamflow and stream temperature implied by the climate change scenarios used in the 5th Assessment Report (AR5) of the IPCC; and (iii) assessed the expected effects of these changes on trout habitat aptitude.To this end, hydrologic simulations with M5 model trees coupled with non-linear water temperature models at the daily time step were fed with high-resolution, downscaled versions of the air temperature and precipitation fields predicted using the most recent climate change scenarios (IPCC, 2013).The effects of basin geology on the stream temperature models and on the estimated changes in thermal regimes were studied.Finally, the changes in the thermal habitat of trout were assessed by studying the violation of the tolerable temperature thresholds of the brown trout.

Materials and methods
The logical framework followed is summarized in Fig. 1.First, the daily global climate models output presented by the IPCC were downscaled to the study area.Then, the obtained local climate models output were applied to generate simulations of streamflow and water temperature.The results are daily values that can be used for the assessment of fish habitat suitability and availability.
The procedure yielded results in the form of continuous time series, but they are presented for two time horizons: the year 2050 (H-2050) and the year 2099 (H-2099).The values for these horizons correspond to the average of the values of the different variables for the decades 2041-2050 and 2090-2099, respectively.

Study sites
In total, 31 sites in 14 mountain rivers and streams inhabited by brown trout were chosen with the aim of encompassing a diverse array of geological and hydrological conditions in the centre of Spain (between the latitudes of 39  1).
The land cover type is mainly pine forest in all of the studied basins (Pinus sylvestris, P. nigra, P. pinea and P. pinaster) (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.Territorial planning does not con-  sider significant changes in land use at mid-century; objectively, changes are not expected after that time because a high percentage of the territory is protected.The studied reaches are not effectively regulated (only small weirs or natural obstacles exist).One large dam lies on the Pirón River (the Torrecaballeros Dam, which has a capacity of 0.32 hm 3 and a maximum depth of 26 m and lies at an altitude of 1390 m a.s.l.), but it does not significantly alter the temporal pattern of streamflow (Santiago et al., 2013).In the Lozoya River, a large dam (the Pinilla Dam, which has a capacity of 38.1 hm 3 and a maximum depth of 30 m and lies at an altitude of 1060 m a.s.l.) exists that separates fish populations above and below the reservoir, although it lies downstream of the studied reach.Hydrological data characterize the streamflow regimes as extreme winter/early spring (groups 13 and 14 in the classification of Haines et al., 1988).However, the hydrographs show a west-to-east smoothing gradient (Fig. 3).This smoothing is associated with the carbonate rocks, whereas greater seasonality is associated with the igneous and detrital geological materials.

Data collection
At each study site, water temperatures were recorded every 2 h throughout the year using 31 Hobo ® Water Temperature Pro v2 (Onset ® ) and Vemco ® Minilog data loggers located at several sites along the studied rivers and streams (Table 1).Loggers were tested for malfunctions before being deployed, and they were placed in areas not exposed to direct sunshine (Stamp et al., 2014).Meteorological data were obtained from nine thermometric and 15 pluviometric sta-Hydrol.Earth Syst.Sci., 21, 4073-4101, 2017 www.hydrol-earth-syst-sci.net/21/4073/2017/ tions of the Spanish Meteorological Agency (AEMET) network, and data from 10 gauging stations (from the official network of the Spanish Water Administration) were obtained to model the streamflows.The AEMET thermometric stations that lie closest to the stream temperature monitoring sites and have at least 30 years of data between 1955 and the present were selected.The selected pluviometric stations were those located within the upstream river basin or near the corresponding gauging station (Table 2).The air temperature and precipitation data from AEMET were tested to assess their reliability by applying a homogeneity test.This test is based on a two-sample Kolmogorov-Smirnov test, and it marks years as possibly containing inhomogeneous data.In the second phase, the marked years are matched against the distribution of the entire series to determine if they contain true inhomogeneities, searching for possible dissimilarities between the empirical distribution functions.Only reliable series were used.The locations of the stations did not change in the studied period.

Climate change modelling and downscaling
Data from nine global climate models associated with the 5th Coupled Model Intercomparison Project (CMIP5) were used, namely BCC-CSM1-1, CanESM2, CNRM-CM5, GFDL-ESM2 M, HADGEM2-CC, MIROC-ESM-CHEM, MPI-ESM-MR, MRI-CGCM3, and NorESM1-M (Santiago et al., 2016).These models provided daily data to simulate future climate changes corresponding to the Representative Concentration Pathways RCP4.5 (a stable scenario) and RCP8.5 (a scenario including a pronounced increase in CO 2 concen-trations) established in Taylor et al. (2009) and used in the AR5 of the IPCC ( 2013).An array of nine general climate models was used to avoid biases due to the particular assumptions and features of each particular model (Kurylyk et al., 2013).Historical simulations of the 20th century were used to control the quality of the procedure and to compare the magnitudes of the predicted changes.Pourmokhtarian et al. (2016) note 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, specifically 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, specifically (1) the speed and (2) direction of the geostrophic wind at 1000 hPa, as well as (3) the speed and (4) direction of the geostrophic wind at 500 hPa.In a second step, the temperature determination was obtained through multiple linear regression analysis using the selected n of the most analogous days.This was performed for the maximum and minimum air temperatures at each station and for each problem day.The linear regression uses forward and backward stepwise selections of the predictors to select only the relevant predictive variables for that particular case.For precipitation, a group of m problem days (the whole days of a month were used) were downscaled together, and the "preliminary precipitation quantity", or the average precipitation of the n most analogous days, was obtained for each prob- lem day.Thus, the m problem days from the highest to the lowest "preliminary precipitation amount" could be sorted.
To assign the final amount of precipitation, each precipitation amount of the m × n analogous days was taken.Then, those m × n amounts of precipitation were sorted, and then those amounts were clustered into m groups.Every quantity was then assigned in order to the m days previously sorted by the "preliminary precipitation amount".Further 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.Such errors are inherently associated with all downscaling methodologies and climate models, which usually introduce bias into their outputs.To eliminate this systematic error, the future climate projections were corrected according to a parametric quantile-quantile method (Monjo et al., 2014), which was performed by comparing the observed and simulated empirical cumulative distribution functions (ECDFs) and linking them using ECDFs obtained from 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, the daily maximum and minimum air temperatures (which were used to infer the mean air temperature) and precipitation were obtained for each climate model, and the whole dataset were used as inputs to simulate the runoff and water temperatures under these climate change scenarios.

Hydrological modelling
Although process-based physical models are considered the standard hydrological models, flexible data-driven machine Hydrol.Earth Syst.Sci., 21, 4073-4101, 2017 www.hydrol-earth-syst-sci.net/21/4073/2017/ learning techniques are gaining popularity because they can be based solely on precipitation and temperature (Shortridge et al., 2016) and can be automatized to perform multiple simulations.Therefore, the prediction of the streamflows under the climate change scenarios was performed with datadriven hydrological models developed using the M5 algorithm (Quinlan, 1992).M5 has been shown to have skill in modelling daily streamflow (Solomatine and Dulal, 2003;Taghi Sattari et al., 2013), including in studies involving climate change (Muñoz-Mas et al., 2016), and it is sufficiently fast to deal proficiently with larger datasets (Quinlan, 2017) Mathematically, M5 is a kind of decision tree that, instead of assigning a single value to each terminal node, assigns a multi-linear regression model (Quinlan, 1992).Therefore, the dataset is hierarchically divided into homogeneous parts and a multi-linear model is adjusted to every part (Hettiarachchi et al., 2005).In this regard, each node, and the corresponding multi-linear regression model, is specialized in particular areas of the data set, such as peak flows or base flows, to name the extremes (i.e. it is a piece-wise linear model with each part dedicated to a particular hydrologic condition) (Taghi Sattari et al., 2013).Based on the multilinear models at the terminal nodes, M5 allows extrapolation, in contrast with other machine learning techniques that have demonstrated little or no extrapolation ability (e.g.random forest or multilayer perceptron) (Hettiarachchi et al., 2005;Shortridge et al., 2016).
The M5 hydrological models were developed in R (R Core Team, 2015) with the Cubist package (Kuhn et al., 2014).One single M5 model tree was trained for each gauging station (ten models were produced in total; Fig. 3 and Table 2), whereas the predictions were supported by the nearest observation (i.e.neighbours = 1) to avoid producing unreliable flows.Finally, M5 was allowed to determine the ultimate number of models (i.e.nodes or areas) into which the dataset is eventually divided (see Kuhn et al., 2014).
Following previous studies, the M5 hydrological models were trained by employing the daily, monthly and quarterly data lags of historical precipitation and air temperature collected at meteorological stations within or nearby the target river basins as input variables (Table 2) (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 and streams that may have different hydrologic behaviours; therefore, the starting set of input variables, which was afterwards subset, was larger than that used in other studies (Solomatine and Dulal, 2003;Taghi Sattari et al., 2013;Muñoz-Mas et al., 2016).The daily variables included the precipitation and air temperature from the current day to the 15th 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 24th previous month (8 variables in total).Consequently, the daily variables overlapped with the current month variable, and the first four quarterly variables overlapped with the monthly data.In the end, 72 variables were gathered, 36 each for air temperature and precipitation.
The whole set of input variables may be relevant for some river systems (Shortridge et al., 2016), although it may cause M5 to overfit the data in others (Schoups et al., 2008).Therefore, the ultimate variable subset was optimized following the forward stepwise approach (Kittler, 1978).This approach relies on iteratively adding input variables (one at a time) while the performance on the test data set improves and stopping (i.e.selecting a smaller subset of the input variables) as soon as the performance stagnates or degrades.However, the classical forward stepwise approach may cause consideration of unrelated variable sets (i.e.disjoint precipitation and air temperature variable lags).To address such potential inconsistencies, the optimization began by testing the precipitation-related variables and only tested the air temperature variables for lags coinciding with those precipitationrelated variables that were already selected.No precautions were taken regarding correlations among inputs (D.P. Solomatine, personal communication, 2016), and the forward stepwise approach sought to maximize the Nash-Sutcliffe efficiency (NSE) index (which ranges from −∞ to 1; Nash and Sutcliffe, 1970) in a fivefold cross-validation (i.e. for each combination of variables, five M5 model trees were trained on four parts and validated with the fifth part, which was held out) (Borra and Di Ciaccio, 2010;Bennett et al., 2013).Finally, in order to account for the uncertainty of the models (Bennett et al., 2013), the variance of the NSEs obtained during the cross-validation was inspected; large intervals led to alternative data partitions.Following previous studies (Fukuda et al., 2013;Muñoz-Mas et al., 2016), once the optimal variables set for each gauging station was determined, 10 M5 model trees (i.e. one per gauging station) were developed using the corresponding subset of variables, and they were used to predict the streamflows under the climate change scenarios.
The daily data were analysed monthly and seasonally using the following statistics: minimum flow (Q min ), the 10th percentile of flow (Q 10 ), the mean flow (Q mean ), and the maximum flow (Q max ).The annual runoff and days of zero flow were also examined.
To assess the significance of the streamflow trends throughout the century, Sen's slope was used (as implemented in the Trend package of R; Pohlert, 2016; pvalue ≤ 0.05) with horizons H-2050 and H-2099.
Finally, the variation in the patterns of the monthly mean streamflow was studied by means of the Ward hierarchical clustering implemented in the cluster R package (Maechler, 2013) on the basis of the rate of change of the normalized monthly mean streamflows in H-2050 and H-2099 and the RCP4.5 and RCP8.5 scenarios.Performance of the obtained cluster was quantified by using the agglomerative coefficient (a.c.).This is a measure of the clustering structure of the dataset, as expressed by Kaufman and Rousseeuw (2005) and its value ranges between 0 (maximum dissimilarity) and 1 (minimum dissimilarity).

Stream temperature modelling
Stream temperature (T s ) at each thermal sampling site was simulated from air temperature (T a ) by means of a modified version of the bounded non-linear regression model described by Mohseni et al. (1998).A previous modification (Term 1 in Eq. 1; Santiago et al., 2016) served to improve the behaviour of the former model, permitting it to be used for daily inputs.In this study, the effect of the instream flow (Q) effect is incorporated.Thus, this model addresses daily mean stream temperature (DMST; T s in Eq. 1) using the daily mean air temperature (DMAT, T a in Eq. 1), the 1-day before variation in the daily mean air temperature ( T a in Eq. 1), and the daily mean flow (Q mean , Q in Eq. 1) as predictors.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 it averages over daily fluctuations in the radiation and heat fluxes.The model is formulated as follows: where µ is the minimum stream temperature ( • C), α is the maximum stream temperature ( • C), β represents the air temperature at which the rate of change of the stream temperature with respect to the air temperature is a maximum ( is the value of the rate of change at β, and λ is a coefficient (dimensionless) that represents the resistance of DMST to change with respect to the 1-day variation in DMAT ( T a ).In the flow component (Term 2 in Eq. 1), ω is the maximum observable variation in stream temperature due to the flow difference (given in • C), τ represents the flow value at which the rate of change of the stream temperature with respect to the flow is a maximum (m 3 s −1 ), and δ (m −3 s) is this maximum rate at τ .Negative values of λ are due to the resistance to stream temperature changes, and thus they must be subtracted from the expected temperature: the more resistant the stream is to temperature change, the closer λ will be to zero.The less resistant the stream is to change, the more negative λ is.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.
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 the Shapiro test and Durbin-Watson test.Moreover, the 7-day lag PACF (partial autocorrelation function) was obtained.These calculations were performed using R. A 95 % confidence interval was calculated for each parameter.Performance was quantified using two indicators: the residual standard error (RSE) and the Nash-Sutcliffe efficiency index (NSE).The Bayesian information criterion (BIC) and the Akaike information criterion (AIC) were used to test the eight-parameter models (Terms 1 and 2 of Eq. 1) against the five-parameter models (Term 1 of Eq. 1).
This model can be classified as semi-physically based.It has some advantages over machine learning methods, such as classification and regression trees (De'ath and Fabricius, 2000) or random forests (Breiman, 2001), because the model parameters imply a mechanistic interpretation of how process drivers act, yielding a higher transferability (Wenger and Olden, 2012).These features make of this model an advantageous option for our goals.

Effects of geology on stream temperature
Geology determines the residence time of deep groundwater in the aquifers underlying streams (Chilton, 1996), and residence times influence discharge temperatures.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 the Bonferroni correction (p-value < 0.05).In the same sense, increments of the annual averages of the daily mean ( T mean ), minimum ( T min ) and maximum ( T max ) stream temperatures were calculated and studied by lithological classes (Table 1).
The variation in the patterns of the monthly mean stream temperature was studied by means of cluster analysis of the temperature increases corresponding to H-2050 and H-2099 for the RCP4.5 and RCP8.5 scenarios (using Ward's hierarchical clustering as implemented in the cluster package of R; Maechler, 2013).As said above, agglomerative coefficient (a.c.) was used as a performance indicator.

Thermal habitat changes
Several tolerance temperatures and thermal niche limits have been described for brown trout (Table 3).The realized niche must reflect energetic efficiency: spending long periods above that threshold makes animals less efficient competitors, and their performance decreases critically (Magnuson et al., 1979;Verberk et al., 2016).Thus, we focused our study on the realized thermal niche.The elected threshold for this study was the occurrence of DMST values above 18.7 • C for seven or more consecutive 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 seven consecutive days is usually the established time for determining thermal tolerance (Elliott and Elliott, 2010), and when this period is exceeded, the death risk (exclusion risk in our case) increases substantially.The chosen threshold was origi- nally determined in one of the streams in this study (the Cega River).
Once DMST was modelled, the frequency of events of seven or more consecutive days above the threshold per year (times above the threshold, TAT ≥ 7), the total days above the threshold per year (DAT), and the maximum consecutive days above the threshold per year (MCDAT) were calculated for the whole period of 2015-2099.
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 current conditions

Longitudinal interpolation and extrapolation
The number of sampling sites and their distribution in the Cega, Pirón and Lozoya rivers (Fig. 2, Table 1) permit the longitudinal interpolation and extrapolation of the predicted water temperatures to study the relationships between the annual average DMST and altitude (strong correlations were detected between these quantities; R 2 = 0.986, 0.985 and 0.881, respectively).A digital elevation model (DEM) with a resolution of 5 m made using lidar and obtained from the National Geographic Institute of the Spanish Government (IGN) was used to perform an altitudinal interpolation of the model parameters to determine the water temperature along the stream continuum to simulate the effects of the climate change scenarios and then to obtain the percentage of stream/river length that will be lost for trout.ArcGIS ® 10.1 software (made by ESRI ® ) was used to manage the DEM.
All variables and abbreviations are summarized in the Appendix A. An overview of the uncertainty issue is given in Appendix B.

Climate change
Under the climate change scenarios, all the meteorological stations will experience noticeable temperature (DMAT) increases through the century.As might be expected, this trend is steeper for the RCP8.5 scenario, especially in summer, though it is also noticeable in winter to a lesser extent (annual trends are shown in Fig. 4; the seasonal results are shown by location in Figs.S1 to S24 in the Supplement).The air temperature variations will run parallel to one another in the two scenarios until mid-century, when the RCP8.5 scenario predicts a similar trend and the increases decrease under the RCP4.5 scenario; the annual change in temperatures for RCP4.5 fluctuates between 2 and 2.5 • C at mid-century and between 2.5 and 3.5 • C at the end of the century (3-4 • C at mid-century and 3.5-4.5 • C at the end of the century in summer) The annual change for RCP8.5 is between 2 and 3 • C at mid-century and between 5 and 7 • C at the end of the century (3.5-4.5 • C at mid-century and 7-8 • C at the end of the century in summer).
The change in the annual precipitation (mm day −1 ) will fluctuate around zero (Fig. 4), although seasonal values will vary (Fig. S1).RCP4.5 predicts a slight decrease (−7 %) by mid-century in total precipitation, which will return to current values by the end of the century.Conversely, RCP8.5

Hydrological regimes
In general, decreases in flow will occur throughout the century, but the degree of change will vary among the sites.Stations located in the western (Tormes) and eastern (Ebrón) extremes of the study area will experience an increase in flow by 2099 after decreasing in the mid-21st century.Lozoya will suffer the most intense flow decreases, followed by Pirón and Cega-Lastras, Tagus and Gallo, and Cabrillas.These patterns of change in flow regimes are predicted to be linked to a westto-east longitudinal gradient; climate change is expected to have less of an influence on discharge at the western stations and Ebrón (in the far eastern portion of the study area).
The hydrological models performed well; all of them achieved NSE values ≥ 0.7 when a number of assorted combinations of variables were selected (Table S1 in the Supplement).Figure 5 shows plots of the monthly Q mean results of the simulations for the RCP4.5 and RCP8.5 scenarios in H-2050 and H-2099.Daily mean streamflow estimated from the climate change model ensemble is given in the Supplement (Dataset S2).

RCP4.5 scenario
Statistically significant (p < 0.05) shifts in the flow regime will be rare in H-2050 (Table 4, Fig. 5).In H-2099, these changes will be less pronounced, but significant changes be-come more frequent (Table 4, Fig. 5).Only two gauging stations (Lozoya and Tagus) exhibit significant reductions in annual discharge.By the end of the century (H-2099), annual discharge is expected to be significantly lower at seven gauging stations.The Tagus Basin will experience the greatest changes in annual discharge.Maximum, mean and minimum daily discharges (Q mean and Q min ), as well as the Q 10 , will become much lower in Tagus River basin.Only Cega-Lastras and Pirón (Duero River basin) will suffer a significant increase in the number of zero-flow days.

RCP8.5 scenario
According to the predictions, the most significant changes in flow regimes will occur at the gauging stations of Cega-Lastras and Lozoya in H-2050 (Table 4, Fig. 5).In H-2099, most sites will experience strong flow reductions, even in seasons where seasonal increases in flow are predicted (e.g.Ebrón and both stations in the Tormes River) (Table 4, Fig. 5).Significant annual runoff reductions in H-2050 will occur at five of the stations, increasing the occurrence of significant losses at 9 out of the 10 sites in H-2099 (i.e.all stations except Ebrón).The most important decreases in every variable and throughout the century were predicted for the stations in the middle Cega Basin and the Tagus Basin.A significant increase in the number of days with no flow was predicted for Cega-Lastras, Pirón and Gallo.

Model parameter behaviour and general trends
The inclusion of the streamflow component improves model performance at 12 out of the 28 study sites (Table 5).In the remaining 16 cases, either no convergence of values was ob- are given in the Supplement (Dataset S3).Performance was high in all the cases except in Pirón 5 where NSE was low.
By the end of the 21st century, the predicted average increase in the mean annual stream temperature among the sites is 1.1 • C for the RCP4.5 scenario (range 0.3-1.6 • C) and 2.7 • C for RCP8.5 (range 0.8-4.0• C).The average increases in maximum annual mean temperature are predicted to be 0.8 • C for RCP4.5 (range 0.1-1.5 • C) and 1.6 • C for RCP8.5 (range 0.2-3.0• C), and the average increases in minimum annual mean temperature are predicted to be 1.0 • C (range 0.4-1.8• C) and 2.7 • C (range 1.1-4.5 • C), respectively.The most important increases are predicted to occur in winter, with summer experiencing smaller increases.

Stream temperature and geological nature
The values of the model parameters showed different behaviours depending on the lithology found in each basin, which thus influences the thermal response to climate change.The thermal amplitude is greater at sites underlain by igneous bedrock (α − µ = 20.38 • C) than at sites underlain by carbonate bedrock (α − µ = 13.07 • C). β values are greater at sites underlain by igneous bedrock (β = 12.71 • C) than at sites underlain by carbonate bedrock (β = 7.80 • C), and λ is significantly greater (λ = −0.140)at sites underlain by carbonate bedrock than at sites underlain by igneous bedrock (λ = −0.292)and at sites underlain by Quaternary detrital material (λ = −0.305)(Fig. 8).All of these differences are significant (p-values < 0.001, as determined using t tests with the Bonferroni correction).
Under the RCP4.5 scenario, T min displays significantly different behaviour at sites underlain by Quaternary detrital material than at sites underlain by carbonate and igneous rocks.Under the RCP8.5 scenario, this difference is solely found between the sites underlain by Quaternary detrital material and those underlain by carbonate rocks.T mean exhibits significant differences between the sites underlain by all three lithologies in both scenarios.All of these results are common to H-2050 and H-2099.In terms of T max , in H-2050, significant differences are found between the sites underlain by carbonate and igneous rocks for both scenarios.These differences are also significant in H-2099 for RCP4.5 and RCP8.5 and between the sites underlain by carbonate rocks and Quaternary detrital material under the RCP8.5 scenario (Fig. 9).
The results of the cluster analysis of the monthly mean stream temperatures revealed a highly homogeneous aggregation of sites for the different combinations of horizons and scenarios, given that the thermal responses of the rivers and streams are tightly linked with lithology (Fig. 7b).The carbonate sites from the Cabrillas stream (in the east) and Pirón 3 (which is strongly influenced by a calcareous spring) form a group of sites that shows low thermal amplitude and in which λ is close to zero.At the other extreme, a group that is made up mainly of sites underlain by igneous material (in the Lozoya and Tormes basins, in addition to several sites found in the detrital basin of Cega-Pirón) shows higher thermal amplitude and lower values of λ than the former group.

The remaining sites have intermediate values of thermal amplitude and resistance.
Table 6.Parameter values of the stream temperature models for every thermograph site (λ is a dimensionless parameter) and values of the performance indicators residual standard error (RSE) and Nash-Sutcliffe efficiency index (NSE).

Effect of streamflow reductions on stream temperature
The predicted flow reductions lead to notable increases in water temperature.The effect of streamflow variation on stream 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 at which the eight-parameter model improves upon the five-parameter model.In all cases, differences in stream temperature between the five-and eight-parameter models are found, and summer flow reductions lead to increases in stream temperature, increasing DAT, TAT ≥ 7 and MCDAT.Among these sites, the threshold is only surpassed at Lozoya and Tormes, increasing the thermal habitat loss.At Cega 1, Cabrillas and Ebrón, α is below the thermal threshold, and at Pirón 1, the stream temperature increase is not sufficient to exceed the threshold.
For all of the sites at which the influence of streamflow on stream temperature was revealed, the eight-parameter model estimates higher values of maximum annual DMST than the five-parameter model.The maximum annual DMST calculated by the eight-parameter model is 3.6 • C higher than that calculated by the five-parameter model at the Tormes 2 site.This difference is not so large at the other sites, and the minimum disagreement between the models (0.01 • C) is noted at the Ebrón and Cabrillas sites.In general, the maximum differences between the two models are noted in igneous catchments, whereas carbonate sites yield the lowest differences.

Effect of climate change on the thermal habitat of brown trout
The length of the thermal habitat of trout will undergo important reductions due to the rises in water temperatures and the increase in the extent of the warm period.In the predictions for H-2050, the 18.7 • C threshold (TAT ≥ 7) will be violated at eight sites under the RCP4.5 scenario and six sites under the RCP8.5 scenario.In H-2099, the threshold will be violated at eight sites under the RCP4.5 scenario and 13 sites under the RCP8.5 scenario.By the end of the century (H-2099), the most notable increases in TAT ≥ 7 (Fig. 10) will be produced at Cega 6, Pirón 5 and Lozoya 3 under the RCP4.5 scenario and at Tormes 1, Cega 4, Cega 6, Lozoya 2, Gallo and Tagus-Poveda under the RCP8.5 scenario.The most significant increases in MCDAT (Fig. 10) will occur at low-altitude sites underlain by igneous rocks and detrital material.In general, the highest temperatures (maximum values of 24.5 • C, Table 7) are predicted to occur in the downstream reaches of the igneous and detrital river basins.In the carbonate basins, only two sites (Tagus-Poveda and Gallo) will exceed the thermal threshold.At mid-century (H-2050), the main changes under the RCP8.5 scenario are similar to those predicted for RCP4.5 at the end of the century (H-2099).RCP4.5 predicts a slower warming from mid-century onwards, whereas RCP8.5 predicts an acceleration of the warming during that period.
Continuous modelling of water temperature by means of the interpolation of model parameters along the Cega, Pirón and Lozoya rivers and the application of the model to DEM data predicts relevant losses of thermal habitat, which will affect up to 56, 11 and 66 % of the lengths of these streams, respectively.In the Cega and Pirón rivers, the habitat loss is expressed relative to the proportion of total stream length where trout currently dwell (98 and 77 km in the Cega and Pirón streams, respectively).In the Lozoya River, the loss is predicted to occur in the reach (20 km) immediately upstream of a large reservoir (the Pinilla reservoir), which produces a total disconnection of the stream.The losses in maximum usable habitat will shift the current downstream limit of the trout distribution from 820 up to 831 m a.s.l in the Pirón River, from 730 up to 830 m a.s.l. in the Cega River, and from 1090 up to 1276 m a.s.l. in the Lozoya River.In the particular case of the Cega River, a window of usable thermal habitat is also predicted to occur upstream from this altitudinal range (from 913 up to 1050 m a.s.l.).

Climate change
Our downscaled results predict greater air temperature increments than the original IPCC (2013) results.These higher temperatures may lead to increased ecological impacts (Magnuson and Destasio, 1997;Angilletta, 2009) caused by the combination of rising water temperatures and decreasing stream flows.The results from the AR5 of the IPCC and its annex, the Atlas of Global and Regional Climate Projections (IPCC, 2013), suggest that droughts are unlikely to increase in the near future for the Mediterranean area.However, air temperatures are expected to rise, subsequently increasing evapotranspiration.As a consequence, the available water in rivers and streams will be reduced.Regional studies have used coarser resolutions than ours, which may be appropriate for their goals (e.g.Thuiller et al., 2006).However, they may be insufficient when more local predictions are needed, as does our study, which treats geographically confined, streamdwelling trout populations.Therefore, fine downscaling techniques like those applied in this study must be used when high-resolution, detailed predictions are needed.

Streamflow
This study predicts significant but diverse streamflow reductions during the present century.At the regional level, a reduction in water resources is expected in 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 of the century (2071-2100), using averages for both the SRES A2 and B1 scenarios (Nakicenovic et al., 2000).Our results predict mean flows that are similar to that value (−23 %; range: 0-49 %), although the emissions scenarios in this study are more severe (that is, they involve greater increases in atmospheric CO 2 ) than those used in the aforementioned studies.
More specifically, the predictions for the RCP4.5 scenario show flow reductions that range from negligibly small to significant (up to 17 %).Under the RCP8.5 scenario, significant reductions become more widespread, ranging up to 49 % of the annual streamflow losses.Our results also predict a relevant increase in the number of days with zero flow for some stations in the detrital area under this scenario (RCP8.5).The predicted streamflow changes are compatible with those obtained in previous studies, although these studies were performed at larger scales (Milly et al., 2005;van Vliet et al., 2013).The apparent differences between the streamflow reductions estimated in this study and those obtained by Milly et al. (2005) and van Vliet et al. ( 2013) (who report lower flow reductions than those given in the present study) might be caused by the regional focus of their predictions (the entire Iberian Peninsula), whereas ours are focused on mountain reaches.
In terms of methods, process-based hydrological models are often preferred for climate change studies (Van Vliet et al., 2012).However, they can be overly complicated and require excessive data inputs, which may also lead to overfitting of the data (Zhuo et al., 2015).Constraining further predictions to within the training domain is a rule of thumb for machine learning studies (Fielding, 1999), although extrapolation is rather common (Elith and Leathwick, 2009).Therefore, taking into account the extrapolation that occurs towards lower flows, which are over-represented in the training dataset, we consider the magnitude of the extrapolation acceptable, and we consider the values, although they are not exempt from uncertainty, to be reliable.

Stream temperature
The model we present in this study showed good performance.Bustillo et al. (2013) recommended the assessment of the impacts of climate change on river temperatures using regression-based methods like ours 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.
However, we also sought to identify relationships between thermal regime and other environmental variables besides 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 temperatures in the remaining 74 % of the cases.Loinaz et al. (2013) quantified the influence of groundwater discharge on temperature variations in the Silver Creek Basin (Idaho, USA), and they concluded that a 10 % reduction in groundwater flow can cause increases of over 0.3 and 1.5 • C in the average and maximum stream temperatures, respectively.Our studied reaches were not influenced by wastewater or reservoir releases (with the exception of releases from the Torrecaballeros Dam on the Pirón River).Kurylyk et al. (2015) showed that the temperature of shallow groundwater influences the thermal regimes of groundwaterdominated 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 models used accurately described the thermal performance of the study sites, and we found significant relationships among the model parameters, the underlying lithologies and the hydrologic responses.Thermal amplitude (α − µ) and temperature at the maximum change rate (β) were lower, and the resistance parameter (λ) was closer to zero, in river basins that were highly influenced by aquifers (mainly carbonate) compared to the others, particularly compared with river basins underlain by carbonate rocks.Since DMST is a variable that is relevant for detecting departures from thermal niche, we can conclude that it is worthwhile to use the more complex eight-parameter model to predict the effects of global warming, especially in igneous catchments.
A wide range of models is described in the literature, and each such model has its strengths and weaknesses.Arismendi et al. (2014) hold that regression models based on air temperature can be inadequate for projecting future stream temperatures because they are only surrogates for air temperature, whereas Piccolroaz et al. (2016) argued that the adequacy depends on the hydrological regime, type of model and the timescale analysis.Their main objections to regressive methods arose when modelling reaches of regulated rivers, but this is not our case.In addition, our model improves the models that were tested in both studies (Arismendi et al., 2014;Piccolroaz et al., 2016).Performance indicators of our models produce good results, showing that the models are sufficiently competent.We show that our model implicitly integrates the effect of other factors, such as geology and flow regime by means of its parameters.A fine mechanistic solution to the modelling issue could need prohibitive methods (Kurylyk et al., 2015), losing the advantages that make attractive the model (input data easy to get).Therefore a compromise between improved precision and increased cost must be met.
The behaviour and dynamics of the parameters offer a promising research field.Their analysis may help to introduce new parametrization criteria to avoid the risk of ignoring the effect of climate warming on groundwater (subsurface water and deep water), for instance.The thermal sensitivity of shallow groundwater differs between short-term (e.g.seasonal) and long-term (e.g.multi-decadal) time horizons, and the relationship between air and water temperatures does not necessarily reflect this difference.This variability should be taken into account in order to avoid un- derestimating the effects of climate warming (Kurylyk et al., 2015).
Regression models are substantially site-specific compared to deterministic approaches (Arismendi et al., 2014).However, the parameters of these regression approaches are still physically meaningful, and these models require fewer variables that can limit the applicability of more complex models 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 describe air temper-atures (and precipitation and streamflow to a lesser extent).On the other hand, our results show that predictions can improve when streamflow is included in the water temperature model, although some streams show little or no sensitivity to the introduction of streamflow into the model.However, the lack of sensitivity is not necessarily due to the absence of the influence of flow on the water temperature but rather to its minor relevance compared to other sources of noise.Thus, when flow data are available, it may be recommended to use the more complex eight-parameter model to predict the effects of climate warming.This conclusion is especially applicable to lithologically sensitive basins, such as those underlain by igneous rocks.
The predicted increase in water temperature will be substantial at most of the study sites.The annual mean rates of change will increase with time.Stewart et al. (2015) predict an increase of 1-2 • C by mid-century in 80 % of the stream lengths in Wisconsin and by 1-3 • C by the latter part of the century in 99 % of the stream lengths, which corresponds to a significant loss in suitable areas for cold-water fish.The results of Stewart et al. (2015) do not differ from ours, except that we expect greater increases by the end of the studied period (up to 4 • C).Our results are also compatible with those of van Vliet et al. (2013), who predicted a water 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 for H-2050 in a river reach in central Spain by midcentury (i.e.daily mean flow reductions between 20 and 29 % and daily mean stream temperature increases up to 0.8 • C).However, we predict that the minima are more sensitive to climate warming than the maxima.

Effects of climate change on brown trout populations
Brown trout are sensitive 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., 2015).In the Iberian Peninsula, the trout distribution is mainly concentrated in mountain streams, where extreme discharges during winter are expected to increase (Rojas et al., 2012).These extreme discharges will likely affect trout recruitment negatively.Thus, the predicted changes in the hydrological regime can subject brown trout populations to more variable conditions, which may occasionally present some populations with insuperable bottlenecks.Trout are polytypic and display an adaptable phenology and rather high intra-population variability in their life history traits that might allow them to show resilience to variations 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 streamflow reduction scenarios is likely to exceed the capacity of many populations to adapt to new conditions (Ayllón et al., 2016).Consistent 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 in the suitable habitat (i.e. the available water volume) (Muñoz-Mas et al., 2016).Finally, the increase in extreme droughts, which involve absolute water depletion, in certain reaches of the streams may be critical for some trout populations.
The predicted increase in winter stream temperatures can affect the sessile phases (i.e.eggs and larvae) of trout development.These phases are very sensitive to temperature changes because it affects their physiology, and because their development is temperature dependent (e.g.Lobón-Cerviá and Mortensen, 2005;Lahnsteiner and Leitner, 2013).Thus, changes in the duration of incubation and yolk sac absorption can affect emergence times and, in turn, the sensitivity of these phases to hydrological regime alterations (Sánchez-Hernández and Nunn, 2016).An increase in stream temperature can also reduce hatchling survival (Elliott and Elliott, 2010).In accordance with the results presented herein, the predicted synergy of streamflow reductions and water temperature increases will cause substantial losses of suitable fish habitat, especially for cold-water fish 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 threshold value) increased by up to 3 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 spatially explicit predictions and to identify reaches where thermal refugia are likely to occur.However, our results suggest that trout will not survive in these reaches because the persistence of thermal refugia is improbable or because their extents will be insufficient.In the Cega, Pirón and Lozoya rivers, important losses of thermal habitat will occur that could jeopardize the viability of the trout population.Behavioural thermoregulatory tactics are 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 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 are able to synergistically contribute to the loss of thermal habitat by increasing daily mean stream temperatures.This effect is especially relevant in summer in the Mediterranean area, when the warmest temperatures and minimum flows usually occur.The existence of thermal refugia represents a possible means of fish survival, and the probability for a water body to become a thermal refugium is highly geologically dependent.In our simulations, the sites that are most dependent on deep aquifers (i.e.basins underlain by Mesozoic carbonate rocks) display improved resistance to 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 south-eastern Spain contain the rear edge of the distribution of native brown trout (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 edge of the trout population in the Iberian Peninsula might shift to the northern mountains to varying extents depending on the presence of relevant mesological features, such as geology.The calcareous mountains of northern Spain could be a refuge for trout because they combine favourable geology and a relatively more humid climate.As a result of this differential response, the western portion of the Iberian range (which is plutonic and less buffered) will eventually experience more frequent local temperature-driven extinction events, thus producing a greater shift northward, than in the eastern Iberian end of this range, which is calcareous and highly buffered and 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 increase less in the central and northern mountains than in the central plateau, and the increases will be smaller in karstic than in granitic (igneous) mountains.At the same time, the side of the peninsula that faces the Mediterranean is expected to be more sensitive to warming and streamflow reductions than the side of the peninsula that faces the Atlantic.Thus, brown trout populations in the karstic mountains of northern Spain (the Cantabrian Mountains and the calcareous parts of the Pyrenees) are better able to resist the climate warming than the populations farther east in the granitic portion of the Pyrenees (Santiago, 2017).Similar patterns may occur in other parts of southern Europe.Most likely, the less pronounced thermal responses of rivers and streams in the karstic areas will allow for greater persistence of the brown trout population, although changes in streamflow regimes will likely also occur there.
In a study of the major basins of Europe, Lassalle and Rochard (2009) predicted that the brown trout would "lose all its suitable basins in the southern part of its distribution area ([the] Black Sea, the Mediterranean, the Iberian Peninsula and the South of France), but [would] likely to continue being abundant in [the] northern basins".Almodóvar et al. (2011) estimated that the brown trout will be eradicated over almost the entire stream length of the studied basins in northern Spain, and Filipe et al. (2013) estimated an expected loss of 57 % of the studied reaches in the Ebro Basin in north-eastern Spain.Our study shows important, yet not so dramatic, reductions in the thermal habitat of Iberian brown trout populations in mountainous areas.The number of general climate models used, the reliability of the downscaling procedure, the resolution of the stream temperature and streamflow models, and the method used to study the threshold imply a substantial improvement in detail (Santiago et al., 2016) over previous work.It is reasonable to infer that many mountain streams appear poised to become refugia for cold-water biodiversity during this century (Isaak et al., 2016).

Conclusions
The main findings of this study are as follows: (i) our downscaled results predict greater air temperature increments than the IPCC's averages, from which our estimations were made; (ii) significant but diverse streamflow reductions are predicted to occur during the present century; (iii) the models presented in this study have been shown to be useful for improving simulations; (iv) the predicted increases in water temperature will be influenced to varying degrees by the flow and geological features of rivers and streams; (v) the thermal habitat of brown trout, a cold-water species, will decrease as a consequence of the synergistic effects of flow reduction and water warming; and (vi) the peaks in water temperature and the complete depletion of the river channels will produce local extinctions, although the ultimate magnitude of the effect will be governed by the geological nature of the basins.
Our findings might be useful in planning the prevention and mitigation of the negative effects of climate change on freshwater fish species at the rear edge of their distributions.A differentiation of areas based on their risk level and viability is necessary to set standardized conservation goals.Our results show that trout conservation requires knowledge 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 provide key information in developing these tools and management options.
In science, and particularly in hydrological studies, the uncertainty is a matter which requires special attention, and it has led us to include a synthesis of our approach to this problem in this appendix.The uncertainty analysis is a necessary step in assessing the risk level in the applicability of a model (Pappenberger and Beven, 2006).Our aim was to study the viability of the brown trout populations, and we modelled the flow and the stream temperature for this purpose.From a conceptual point of view our approaches were consistent with it.
On data inputs to build the models, uncertainties and inconsistencies are a habitual issue (Juston et al., 2013).The meteorological and hydrological services subject data to their own quality controls but systematic error cannot always be completely controlled (Beven and Westerberg, 2011;McMillan et al., 2012).For this reason, in addition, we tested the input data seeking inconsistencies.
The modelling of the river reaches as one-dimensional elements implies a simplification of the fluvial ecosystem that is generally accepted at this scale (e.g.Viganò et al., 2015;Ahmed and Tsanis, 2016), especially for ecological purposes (e.g.Caiola et al., 2014).Nevertheless, the size of the rivers under study made little or nothing relevant the variations in width and depth (it was verified in the field).
Regarding the parameterization of the models, crossvalidation was used to evaluate the uncertainty in these process, and indicators such as the NSE (for hydrological and thermal models), the deviance and the RSE (for thermal models) were calculated.In the case of the thermal model, the functions of distribution of the parameters of the model were built by non-parametric bootstrap, and the mean values were chosen as the most proficient estimators.As results show, pa-rameters tell us about the functional behaviour of catchments (particularly on the effects of the catchments geology on the streams temperature) and this might improve predictions in ungauged basins by better controlling uncertainty (Juston et al., 2013).
Once the models were constructed, It was verified that the overlaps of the ranges of the model input variables and the ranges of the outputs were significant (p < 0.05).The non-overlapping zones affected, on the one hand, infrequent events (great floods) and the extreme temperature zone (zone of extrapolation), being the last one the scope in which we expected to work.However, the weakness of the hydrological model in the flood zone should be considered for other applications and developments of the model.The hydrological model given us sufficient and relevant information since its possible weak points (extrapolation in the floods assessment) did not affect our goal: we focused on central trends and minimum values, and they were solidly represented in the samples and in the simulations.
As said, the inherent uncertainty of the climate predictions according to the scenarios RCP4.5 and RCP8.5 was attenuated by means of the ensemble technique, showing the dispersion of the results by mean of the percentiles in Figs.S1 to S24 (in the Supplement).Beven (2011) exposed his legitimate concerns on the credibility of climate models which fail when are compared with the control period and, consequently, we used ERA-40 reanalysis to control this source of bias with excellent results.
There was a time dependence between the errors of the model and the scope of the prediction, but these errors were only important in the zone of high temperature and low flow, as expected by the physical nature of the climatic variables.Moreover, this is the behaviour of the variables that was our intention to evaluate.

Figure 1 .
Figure 1.Logical framework of the study.

Figure 2 .
Figure 2. River network and location of the study sites (water temperature data loggers), with details regarding lithology.The grid depicts the actual occurrence of brown trout in Spain.
• 53 and 41 • 21 N).Specifically, the investigated sites are located in the Tormes River and its tributaries, the Barbellido River, the Gredos Gorge and the Aravalle River (in the Duero Basin); the Cega River and the Pirón River (the Pirón River is a tributary of the Cega River in the larger Duero Basin); the Lozoya River, the Tagus River, the Gallo River, and the Cabrillas River (all four of which are in the Tagus Basin); the Ebrón River and the Vallanca River (the Vallanca River is a tributary of the Ebrón River in the Turia Basin); and the Palancia River and the Villahermosa River (Fig.2).The ma-jor geological components that lithologically characterize the mountain sites in the Duero and Lozoya basins are igneous rocks; the altitudinally lower sites in the Duero Basin are underlain by Cenozoic detrital material, and the eastern basins (Tagus, Gallo, Cabrillas, Ebrón, Vallanca, Palancia and Villahermosa) are underlain primarily by Mesozoic carbonates.The distribution of geological materials was retrieved from the Lithological Map of Spain (IGME, 2015) (Table

Figure 3 .
Figure 3. River regime patterns for the different gauging stations.The flows are expressed as percentages of the mean annual flow, and the months (horizontal axis) are ordered from January to December.

Figure 4 .
Figure 4. Changes in mean air temperature and total annual precipitation related to climate change for the nine general climate models and the two climate change scenarios for the all the studied meteorological stations.

Figure 5 .
Figure 5. Predicted monthly mean specific flow in H-2050 and H-2099 for the RCP4.5 and RCP8.5 scenarios.Shaded areas indicate decadal fluctuations.Triangles show significant negative or positive trends (Sen's slope p ≤ 0.05); the sign of each trend is indicated by the directions in which the triangles point.

Figure 7 .
Figure 7. Study sites clustered by the predicted change ratios of the seasonal mean streamflow (gauging stations) and by the predicted increase in the monthly mean stream temperature ( • C) at the water temperature recording sites in H-2050 and H-2099 for the RCP4.5 and RCP8.5 scenarios.Axes indicate geographic positions (UTM coordinates).The colours and numbers indicate the clusters.

Figure 8 .
Figure 8. Distributions of the stream temperature model parameter values (α − µ, β, γ and λ) in relation to lithology.Differences were assessed using Student's t test with the Bonferroni correction (p < 0.05).

Figure 9 .
Figure 9. Distributions of T min , T mean and T max in relation to lithology for the climate change scenarios RCP4.5 and RCP8.5 in H-2050 and H-2050.The reference period corresponds to the simulated period 2010-2019.

Figure 10 .
Figure 10.Increases in TAT ≥ 7 (time above the threshold during seven or more consecutive days), MCDAT (maximum consecutive days above the threshold) and DAT (days above the threshold per annum) from the present to H-2050 and H-2099 for RCP4.5 and RCP8.5.

Table 1 .
Description of the data logger (thermograph) sites, specifying given name, UTM coordinates (Europe WGS89), altitude (metres above 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 (lithology) of the data logger site (the last of which was obtained from IGME, 2015).Bold letters indicate sites associated with the gauging stations.

Table 2 .
Official stations used (meteorological and hydrological), variables, length of time series used and geographical position.AEMET: Spanish Meteorological Agency; CHD: Water Administration of Duero Basin; CHT: Water Administration of Tagus Basin; and CHJ: Water Administration of Júcar Basin.

Table 3 .
Different classes of thermal thresholds for emerged trout classes found in the literature.The type of experiment differentiates the experiments with controlled (laboratory) and uncontrolled (wild) temperature.Latitude of the experiments' location is shown.

Table 5 .
Bayesian (BIC)and Akaike (AIC) information criteria values for the stream-temperature models with five and eight parameters.

Table 7 .
Maximum daily mean stream temperature ( • C) in each site in the year 2015 and the horizons H-2050 and H-2099.Both scenarios (RCP4.5 and RCP8.5) are shown.Bold numbers: values > 18.7 • C.