Climate Trends Impact on the Snowfall Regime in Mediterranean Mountain Areas: Future Scenario Assessment in Sierra Nevada (Spain)

Snow constitutes a key component of the water cycle, which is directly affected by changes in climate. Mountainous regions, especially those located in semiarid environments, are highly vulnerable to shifts from snowfall to rainfall. This study evaluates the influence of future climate scenarios on the snowfall regime in the Sierra Nevada Mountains, an Alpine/Mediterranean climate region in southern Spain. Precipitation and temperature projections from two future climate scenarios representative concentration pathway (RCP) 4.5 and RCP 8.5, Fifth Assessment Report of the Intergovernmental Panel for Climate Change (AR5 IPCC)) were used to estimate the projected evolution of the snowfall regime on both annual and decadal scales during the period of 2006–2100. Specific snowfall descriptors of torrentiality are also analyzed. A general decrease of the annual snowfall was estimated, with a significant trend that ranged from 0.21 to 0.55 (mm·year−1)·year−1. These changes are dependent on the scenario and region in the study area. However, the major impact of future climate scenarios on the snowfall regime relates to an increased torrentiality of snowfall occurrence, with a decreased trend of the annual number of snowfall days (RCP 4.5: −0.068 (days·year−1)·year−1 and RCP 8.5: −0.111 (days·year−1)·year−1) and an increased trend in the annual mean snowfall intensity (RCP 4.5: 0.006 (mm·days−1)·year−1 and RCP8.5: 0.01 (mm·days−1)·year−1)) under both scenarios. This enhanced torrentiality is heterogeneously distributed, with the most semiarid region, which is currently the one least influenced by snow, being the region most affected within the study area.


Introduction
Snow is a basic component of the Earth's surface energy balance [1]. Climate trends, especially short-term and long-term precipitation and temperature variations, condition the snow regime and the hydrological evolution of the river flow response at the basin and sub-basin scales. The link between both climate and flow trends is crucial in mountainous areas, where small variations in temperature, among others, can produce significant impacts on precipitation (which occur as either rainfall or snowfall), snowmelt, and sublimation/evaporation. This can consequently result in very different flow signatures [2][3][4].
In recent years, different authors have concluded that rainfall patterns are changing all over the globe as a result of global warming [5][6][7][8]. Moreover, the analysis of precipitation events shows an increase in their spatiotemporal variability during the last decades [9][10][11]. These variations, together with the increase in temperature, affect both the amount and distribution of snowfall, especially in semiarid regions [12]. The impacts of this climate variability on the snow and hydrological regime are more evident over semiarid high mountain regions (i.e., areas higher than 1500 m a.s.l.) due their particularly extreme conditions [13,14], in which the high variability of the annual and seasonal climate regimes is usually propagated to and amplified by the river flow. This is the case in the Mediterranean mountain areas, where alpine and semiarid conditions coexist [15].
The impact of climatic variability in the medium and long term (i.e., 20-and 100-year time horizons, respectively) hydrological regimes in such areas is being assessed from different approaches on different spatiotemporal scales [14,16,17]. The development of local climate indicators from both climate and hydrological variables are key to addressing the economic impacts of the expected changes and to assessing decision-making processes [18][19][20][21].
This work assesses the influence of future climate scenarios on the snow regime in the Sierra Nevada Mountains, an alpine climate region in South Spain, during the 21st century. The Sierra Nevada area is a national park and a biosphere reservoir, with altitudes ranging from 2000 to 3500 m a.s.l. The park is a representative example of snow regions in semiarid environments and is characterized by a high inter-and intra-annual climate variability that influences most of the hydrological processes on different time and spatial scales [22,23]. For this reason, Sierra Nevada (Spain) is part of the global climate change observatories network [24]. The future evolution of climate and its impact on snow occurrence and persistence not only affects the hydrological regime [23,25] in the area but also poses a risk for the ecosystems associated with the snow domain and the numerous endemic species identified in the park [26,27].
For this assessment, precipitation and temperature projections in this area from two future climate scenarios derived from a Global Circulation Model (CGM) were used to estimate the evolution of the snowfall regime on both annual and decadal scales in the 2006-2100 period. Additionally, specific torrentiality indicators, which are defined in Section 2, were analyzed together with the spatial distribution of the resulting impacts.

Study Area
Sierra Nevada ( Figure 1) is a mountainous range located in southeastern Spain (37 • 6'0" N; 3 • 5'24" W), approximately 80 km long in the east-west direction and 27 km wide in the north-south direction. It runs parallel to the Mediterranean Sea coastline, approximately 40 km from its highest peaks (3479 m a.s.l.). The area has a characteristic steep topography, with strong altitudinal gradients and a large variety of associated ecosystems. The annual precipitation regime is highly variable and ranges between values close to 1000 mm in wet years and 200 mm in dry years [12]. Snow recurrently appears in the mountain range in altitudes higher than 1000 m a.s.l. and is mostly persistent over 2000 m a.s.l. from November to May, undergoing several cycles of accumulation-ablation during the snow season. This results in a very heterogeneous spatial distribution over the years. The mean annual fractional snow cover area for the period 2000-2013 was 0.21 m 2 ·m −2 , ranging from 0.9 to 0.16 m 2 ·m −2 in wet/cold and dry/warm years, respectively, with a mean standard deviation of 0.23 m 2 ·m −2 [28]. This behavior makes Sierra Nevada a relevant biodiversity enclave [26,29,30].
For this study, the area is divided into five regions belonging to the headwaters of the main river basins in Sierra Nevada: Adra (R1), Andarax (R2), Fardes (R3), Genil (R4), and Guadalfeo (R5). These regions are representative of the N-S/W-E gradients in the area, with R5, R1, and R2 being influenced by the sea and the decreasing precipitation regime during the last two decades. R4 and R3 correspond, respectively, to the most continental and semiarid regions in the area, which have a more stable precipitation regime. Table 1 shows some relevant climate descriptors of these regions during the reference period   [12,31], as well as the associated regional annual trends. The weather station network used in this study for measuring precipitation (yellow circle) and temperature (yellow square), and the limits of the five regions into which the study area was divided for the spatial analysis: Adra (R1), Andarax (R2), Fardes (R3), Genil (R4), and Guadalfeo (R5). Table 1. Summary of the selected characteristics of each region (see Figure 1) in the study area (SN), together with the observed mean annual trends for the climate variables during the reference period 1961-2000. The trend analysis results are defined according to the Mann-Kendall test as having either (1) excellent significance, 99% (***); (2) good significance, 95% (**); (3) acceptable significance, >90% (*); or (4) no significance, less than 90% (+) [12,31].  The weather station network used in this study for measuring precipitation (yellow circle) and temperature (yellow square), and the limits of the five regions into which the study area was divided for the spatial analysis: Adra (R1), Andarax (R2), Fardes (R3), Genil (R4), and Guadalfeo (R5). Table 1. Summary of the selected characteristics of each region (see Figure 1) in the study area (SN), together with the observed mean annual trends for the climate variables during the reference period 1961-2000. The trend analysis results are defined according to the Mann-Kendall test as having either (1) excellent significance, 99% (***); (2) good significance, 95% (**); (3) acceptable significance, >90% (*); or (4) no significance, less than 90% (+) [12,31].  Figure 2 shows the scheme followed in assessing the impacts of future climate scenarios on the snow regime. Information was available regarding the reference values of the climate variables analyzed in this study during the reference period 1961-2000 from previous works [12,31], on both annual and decadal scales. T: temperature, P: total precipitation (rainfall and snowfall), S: snowfall, E: number of snowfall events per year, and D: annual mean duration of snowfall events. Figure 2 shows the scheme followed in assessing the impacts of future climate scenarios on the snow regime. Information was available regarding the reference values of the climate variables analyzed in this study during the reference period 1961-2000 from previous works [12,31], on both annual and decadal scales. AEMET is the Spanish Meteorological Agency, which is the data source of the historical and projected climate variables for future scenarios RCPs4.5 and 8.5 in the Fifth Assessment Report of the Intergovernmental Panel for Climate Change (AR5 IPCC)) [14]. CIIs are the climate impact indicators defined for torrentiality assessment.

Data Sources
Precipitation and temperature data series were provided by the Spanish Meteorological Agency (AEMET) in specific locations across the study area for the reference  and future (2006-2100) periods. Four different representative concentration pathways (RCPs) were available in the framework of the Fifth Assessment Report of the Intergovernmental Panel for Climate Change [14], which assesses different levels of greenhouse gas emissions in future years. For this study, two of these scenarios were chosen: an intermediate mitigation scenario (RCP 4.5), which assumes an emission reduction from 2050 onwards, and a high emission trajectory (RCP 8.5), which contemplates the possibility of not reducing emissions.
Datasets obtained from different global circulation models (GCM4, MPI-ESM-LR, ECHO, GFDL, GISS-AOM, GISS-ER, INM-CM4, MIR-GCM2-0, and MIR-GCM2-1) were available from two downscaling techniques each (analog and statistical downscaling methods system (SDMS)). AEMET is the Spanish Meteorological Agency, which is the data source of the historical and projected climate variables for future scenarios RCPs4.5 and 8.5 in the Fifth Assessment Report of the Intergovernmental Panel for Climate Change (AR5 IPCC)) [14]. CIIs are the climate impact indicators defined for torrentiality assessment.

Data Sources
Precipitation and temperature data series were provided by the Spanish Meteorological Agency (AEMET) in specific locations across the study area for the reference  and future (2006-2100) periods. Four different representative concentration pathways (RCPs) were available in the framework of the Fifth Assessment Report of the Intergovernmental Panel for Climate Change [14], which assesses different levels of greenhouse gas emissions in future years. For this study, two of these scenarios were chosen: an intermediate mitigation scenario (RCP 4.5), which assumes an emission reduction from 2050 onwards, and a high emission trajectory (RCP 8.5), which contemplates the possibility of not reducing emissions.
Datasets obtained from different global circulation models (GCM4, MPI-ESM-LR, ECHO, GFDL, GISS-AOM, GISS-ER, INM-CM4, MIR-GCM2-0, and MIR-GCM2-1) were available from two downscaling techniques each (analog and statistical downscaling methods system (SDMS)). Although the use of the multi-model ensemble mean (MMEM) of climate model runs is generalized in climate change impact studies, some authors have pointed out that there is no reason to expect that the MMEM always provides the best estimate of climate change signals [32]. Even if it is somewhat obvious that an ensemble of models should provide more information than a single one, it is not necessarily clear how to optimally use that information in practice [33]. Different authors [32,34,35] have proposed the use of metrics to choose among climate models in performing impact studies when the use of the MMEM is not likely to improve future climate assessments. In this context, metrics refer to quantitative indicators of model performance during past or present-day climate simulations against the available observations, under the assumption that this metric is expected to also be informative on the likelihood that the model will correctly simulate the target variable under a future climate scenario.
Taking into account the extension and steep topography of the study area and the spatial resolution of climate models, a first assessment of the available downscaled modeled projections was conducted prior to the impact analysis. Two metrics were defined after the analysis of the historical data series at the local scale: (1) the annual mean of the daily minimum temperature (Tmin) and (2) the annual precipitation (P). Each combination of model and downscaling was evaluated in each weather station in the area during the reference period, and the root squared mean error (RMSE) was computed. The analog method resulted in lower RMSE values for Tmin in all cases. However, in only two from the whole set of models, these values were lower than one degree (one of them was close to this value, however), which is a constraint for snowfall estimation. Regarding P, similar RMSE values were found between analog and statistical methods of downscaling, but differences were observed up to 50% among models. Focusing on the two models selected by Tmin RMSE, the difference of the P RMSE value was close to 20%. From these results, it was clear that among the 18 different combinations of modeling and downscaling, the Max Planck Institute Earth system model in low resolution (MPI-ESM-LR) [36] downscaled using the analog method [37] was the combination that best performed in the study area (RMSE = 0.8 C and RMSE = 201 mm for Tmin and P, respectively). Thus, this combination was selected for the present work. The global model selected, MPI-ESM-LR, is a coupling between the ECHAM6 [38] for the atmosphere and the MPI Ocean Model (MPI-OM) [39] for the oceans [40][41][42][43] at low resolution. The selected analog downscaling method regionalizes the GCM results based on the use of synoptic analogs for each day in the time series [44,45]. From this combination, daily data were obtained from a total of 44 precipitation and four temperature weather stations in the Sierra Nevada area ( Figure 1).

Interpolation of Meteorological Variables
Specific spatial interpolation algorithms that include the effects of topography [46] were used to distribute the point meteorological daily information throughout the study area [47].
For each variable (i.e., precipitation and temperature) a linear gradient with elevation was locally defined from the historical datasets using a least square error method. Afterwards, the residuals between each station value that were obtained using the linear gradient were calculated. These residuals were interpolated by square inverse distance weighting (IDW2) by using the three closest stations [48,49] and added to the theoretical value calculated using the linear gradient [50,51]. This procedure had been previously used in the Sierra Nevada area for other climate studies [12,[52][53][54]. The spatial resolution for the gridded resulting map series is 30 m.

Snowfall Analysis and Snow Climate Indicators
Although the constant threshold temperature approach is a simple method, with snowfall occurrence being a more complex process depending on several atmospheric variables [55] and the density and size of snow particles [56], the available dataset did not allow a different approach to be performed, as other works have addressed [57][58][59]. However, the fixed threshold temperature approach has been and is still used in many works [60][61][62], including in the study area [54,63], and it has proven to be appropriate on a daily scale, especially for hydrological applications [64].
Snowfall occurrence and amount was determined on a daily basis from the use of a threshold temperature, T SN , over which all precipitation is considered rainfall, following Equation (1). In this equation, "Snowfall" is the daily amount of snowfall, P the daily amount of total precipitation occurring below T SN , and T is the temperature during the day. For each day, the temperature datasets include maximum, minimum, and mean temperature values. A triangular shape was assumed for sub-daily variations and the fraction of time with T below T SN was adopted as the snowfall fraction of the daily precipitation. In this work, a constant value of 0.5 • C was assumed for T SN , which has been proven to perform satisfactorily in previous studies [28,53,63].
Water 2018, 10, 720 6 of 22 From the resulting daily series of snowfall following this approach, two climate impacts on snow indicators (CIIs) were used to quantify the mean intensity of snowfall events and their torrentiality: • Snowfall days indicator (Dsnowfall): number of days each year in which snowfall is higher than 0.5 mm. • Snowfall intensity (Isnowfall): relationship between the annual snowfall and the snowfall days indicator, as outlined by Equation (2).

Long-Term Data Analysis
From the resulting projected 30 × 30 m maps of the downscaled daily weather and snowfall variables in Sierra Nevada during the 2006-2100 period, the analysis on the annual and decadal scales was carried out by aggregation of the daily data associated with the selected scenarios, RCP 4.5 and RCP 8.5. The results were also averaged over both the whole study area and each of the five regions identified in this work. These results were annual precipitation (P); annual mean values of the mean, maximum, and minimum daily temperature (Tmean, Tmax, Tmin); and annual snowfall (Snowfall). Additionally, the defined CIIs were also assessed for each region and time scale from the snowfall map time series, Dsnowfall and Isnowfall, in order to identify the regions which will likely be the most influenced by the estimated shifts of the snowfall regime.
Trend analysis was performed on both the resulting annual and decadal time series of the different variables for the future time study period, which was expressed as absolute anomalies from the respective mean value during the reference period (Table 1), with significance assessment following the Mann-Kendall test [65]. Autocorrelation testing [66] was previously performed to guarantee consistency in the results. Figure 3 represents the annual evolution of the average daily mean, maximum, and minimum temperature in the Sierra Nevada area for the future projections during 2006-2100, expressed as absolute anomalies from the mean value during the reference period (Table 1). On an annual basis, an increasing trend can be observed in all variables for both scenarios analyzed. RCP 8.5 was the most adverse scenario in terms of warming, as expected. The mean temperature anomaly follows a slightly variable regime during the future period in both scenarios, being positive in all cases. The maximum and minimum temperature anomalies exhibit a highly variable regime, with negative values in both scenarios. This can also be observed in the decadal regime. The resulting trends are significant on a high level of confidence (99%) in all cases. Autocorrelation was found only in the case of the annual mean values of the daily minimum and mean temperature under RCP 8.5 scenario. Therefore, the resulting significance in the trend of these variables can be affected by this fact.

Long-Term Precipitation and Temperature Analysis for Future Scenarios
Negative anomalies for both the annual and decadal mean values of the daily maximum and minimum temperature can be found under the RCP 4.5 scenario, and they expand over the whole future study period, with some apparent stabilization of the decadal mean daily maximum and mean temperature during the last part of the period. This is not observed under the severe scenario of RCP 8.5, with anomalies close to 6 • in both variables, on both the annual and decadal scales, at the end of the study period.

fact.
Negative anomalies for both the annual and decadal mean values of the daily maximum and minimum temperature can be found under the RCP 4.5 scenario, and they expand over the whole future study period, with some apparent stabilization of the decadal mean daily maximum and mean temperature during the last part of the period. This is not observed under the severe scenario of RCP 8.5, with anomalies close to 6° in both variables, on both the annual and decadal scales, at the end of the study period.  The annual and decadal mean daily minimum temperature anomalies show the highest variability in all cases. Under the RCP 8.5 scenario, this variable follows a trend towards prevailing positive values during the last decades of the study period. Figure 4 shows the annual and decadal evolution of the average annual precipitation in the Sierra Nevada area for future projections during the study period from 2006 to 2100 expressed as absolute anomalies from the mean value during the reference period ( Table 1). As opposed to temperature trends, a decreasing evolution is observed in both scenarios. However, on a significant level, an annual trend is found only for the RCP 8.5 scenario, the decreased trend in precipitation of which is four-fold that found in the RCP 4.5 scenario. It must be noted that the mean annual precipitation in the period differs 4% between both scenarios, being 503 and 482 mm·year −1 for RCP Water 2018, 10, 720 8 of 22 4.5 and RCP 8.5, respectively. However, a large difference is found for the associated standard deviation (i.e., 44.11 and 47.37 mm·year −1 for RCPs 4.5 and 8.5, respectively).
positive values during the last decades of the study period. Figure 4 shows the annual and decadal evolution of the average annual precipitation in the Sierra Nevada area for future projections during the study period from 2006 to 2100 expressed as absolute anomalies from the mean value during the reference period ( Table 1). As opposed to temperature trends, a decreasing evolution is observed in both scenarios. However, on a significant level, an annual trend is found only for the RCP 8.5 scenario, the decreased trend in precipitation of which is four-fold that found in the RCP 4.5 scenario. It must be noted that the mean annual precipitation in the period differs 4% between both scenarios, being 503 and 482 mm·year −1 for RCP 4.5 and RCP 8.5, respectively. However, a large difference is found for the associated standard deviation (i.e., 44.11 and 47.37 mm·year -1 for RCPs 4.5 and 8.5, respectively).  Decadal trends are not significant, however. In fact, the first three decades of RCP 8.5 and the first four of RCP 4.5 reflect an increasing evolution, which results in a non-significant trend (i.e., no trend) in RCP 4.5. Decadal anomalies are larger and more fluctuating in the case of RCP 8.5.
These results have a clear impact on the snowfall occurrence distribution, as described in the next section. This is particularly relevant in the highest areas, as Figure 5 shows, where both Decadal trends are not significant, however. In fact, the first three decades of RCP 8.5 and the first four of RCP 4.5 reflect an increasing evolution, which results in a non-significant trend (i.e., no trend) in RCP 4.5. Decadal anomalies are larger and more fluctuating in the case of RCP 8.5.
These results have a clear impact on the snowfall occurrence distribution, as described in the next section. This is particularly relevant in the highest areas, as Figure 5 shows, where both increased temperature and decreased precipitation mean values for the whole study period can be observed. increased temperature and decreased precipitation mean values for the whole study period can be observed.

Long-Term Snowfall Analysis for Future Scenarios
Based on the previous results, daily snowfall occurrence was estimated for each scenario on a 30 × 30 m grid, as described in Section 2. In this work, the results were aggregated on the annual and decadal scales and averaged over the different spatial units in this work. Figure 6 shows the spatial distribution of the anomalies of the mean annual snowfall in the study area during the future period under both climate scenarios in this work, together with the mean annual trends during the study period. The largest impacts can be found in the highest areas (i.e., mostly negative anomalies), being both larger and more extended in RCP 8.5, as expected. Only under RCP 4.5 did some areas in the west side of the Sierra Nevada basins exhibit positive anomalies, although these were much lower in absolute values than the negative ones. Trends are negative in all cases, and they are higher in the snow domain area.

Long-Term Snowfall Analysis for Future Scenarios
Based on the previous results, daily snowfall occurrence was estimated for each scenario on a 30 × 30 m grid, as described in Section 2. In this work, the results were aggregated on the annual and decadal scales and averaged over the different spatial units in this work. Figure 6 shows the spatial distribution of the anomalies of the mean annual snowfall in the study area during the future period under both climate scenarios in this work, together with the mean annual trends during the study period. The largest impacts can be found in the highest areas (i.e., mostly negative anomalies), being both larger and more extended in RCP 8.5, as expected. Only under RCP 4.5 did some areas in the west side of the Sierra Nevada basins exhibit positive anomalies, although these were much lower in absolute values than the negative ones. Trends are negative in all cases, and they are higher in the snow domain area.
Water 2018, 10, x FOR PEER REVIEW 9 of 22 increased temperature and decreased precipitation mean values for the whole study period can be observed.

Long-Term Snowfall Analysis for Future Scenarios
Based on the previous results, daily snowfall occurrence was estimated for each scenario on a 30 × 30 m grid, as described in Section 2. In this work, the results were aggregated on the annual and decadal scales and averaged over the different spatial units in this work. Figure 6 shows the spatial distribution of the anomalies of the mean annual snowfall in the study area during the future period under both climate scenarios in this work, together with the mean annual trends during the study period. The largest impacts can be found in the highest areas (i.e., mostly negative anomalies), being both larger and more extended in RCP 8.5, as expected. Only under RCP 4.5 did some areas in the west side of the Sierra Nevada basins exhibit positive anomalies, although these were much lower in absolute values than the negative ones. Trends are negative in all cases, and they are higher in the snow domain area.   Table 2 shows the mean values of the variables under analysis during the study period, averaged over each region in Sierra Nevada and the whole study area (SN) under both future climate scenarios, together with the results of the trend analysis performed for the whole study period. Table 2. Summary of the mean values of the study variables during the 2006-2100 study period under each scenario, averaged over each region (see Figure 1) in Sierra Nevada (R1 to R5) and the whole study area (SN), together with the associated mean annual trends. The trend analysis results are defined according to the Mann-Kendall test as having either (1) excellent significance, 99% (***); (2) good significance, 95% (**); (3) acceptable significance, >90% (*); and (4) no significance, less than 90% (+).    The decadal evolution of the mean annual snowfall for the future period under each scenario can be observed in Figure 8 as anomalies, averaged for each region in Sierra Nevada and the whole study area. The decadal regime more efficiently shows the decrease of snowfall associated with both scenarios in all regions, with the RCP 8.5 scenario (being the most severe) causing the highest impact in snowfall occurrence in all cases. Annual variability is large in all projected decades, with the anomalies in R2 being lower than those in the rest of the regions. However, this region is also the most semiarid area in Sierra Nevada, with an average mean annual snowfall of 44 mm·year −1 during the 1961−2000 reference period, which makes the impact of future climate scenarios highly relevant for all ecotypes associated with snow occurrence in this region.
Anomalies are positive in most regions for the first half of the future period under the RCP 4.5 scenario; however, they drop to negative values in all regions after three decades under the RCP 8.5 scenario. Figure 9 shows the annual evolution of the Dsnowfall indicator defined in Section 2.5 during the 2006-2100 period for each region in Sierra Nevada (Figure1) and the whole study area. This variable is also related to the annual mean intensity of snowfall (Isnowfall) that is shown in Figure  10. Again, it is illustrated that R5, with higher Dsnowfall values (ranging between 44 and 8 days/year under the RCP 4.5 scenario and between 38 and 4 days/year under the RCP 8.5 scenario), is the most snow-influenced domain in Sierra Nevada. The semiarid character of R2 is also highlighted (ranging between 26 and 1 days/year under the RCP 4.5 scenario and 18 and 0 days/year under the RCP 8.5 scenario). Significant decreasing trends (99%) are found for both scenarios in all cases; however, relevant differences are found between scenarios.
Regarding Isnowfall (Figure10), this indicator shows similar values in most of the identified regions in Sierra Nevada, with the exception of R2 and R3, which exhibit more torrential behavior. The decadal evolution of the mean annual snowfall for the future period under each scenario can be observed in Figure 8 as anomalies, averaged for each region in Sierra Nevada and the whole study area. The decadal regime more efficiently shows the decrease of snowfall associated with both scenarios in all regions, with the RCP 8.5 scenario (being the most severe) causing the highest impact in snowfall occurrence in all cases. Annual variability is large in all projected decades, with the anomalies in R2 being lower than those in the rest of the regions. However, this region is also the most semiarid area in Sierra Nevada, with an average mean annual snowfall of 44 mm·year −1 during the 1961−2000 reference period, which makes the impact of future climate scenarios highly relevant for all ecotypes associated with snow occurrence in this region.
Anomalies are positive in most regions for the first half of the future period under the RCP 4.5 scenario; however, they drop to negative values in all regions after three decades under the RCP 8.5 scenario. Figure 9 shows the annual evolution of the Dsnowfall indicator defined in Section 2.5 during the 2006-2100 period for each region in Sierra Nevada (Figure1) and the whole study area. This variable is also related to the annual mean intensity of snowfall (Isnowfall) that is shown in Figure 10. Again, it is illustrated that R5, with higher Dsnowfall values (ranging between 44 and 8 days/year under the RCP 4.5 scenario and between 38 and 4 days/year under the RCP 8.5 scenario), is the most snow-influenced domain in Sierra Nevada. The semiarid character of R2 is also highlighted (ranging between 26 and 1 days/year under the RCP 4.5 scenario and 18 and 0 days/year under the RCP 8.5 scenario). Significant decreasing trends (99%) are found for both scenarios in all cases; however, relevant differences are found between scenarios.
Regarding Isnowfall (Figure10), this indicator shows similar values in most of the identified regions in Sierra Nevada, with the exception of R2 and R3, which exhibit more torrential behavior. As previously stated, these regions represent the semiarid conditions in Sierra Nevada, and their features tend to be enhanced in the future climate scenarios, especially under RCP 8.5, as expected. All resulting trends are positive (RCP 4.5: 0.006 (mm·days −1 )·year −1 and RCP 8.5: 0.01 (mm·days −1 )·year −1 ), although quite moderate under RCP 4.5. The decadal scale more clearly confirms the semiarid character of R2, with high variability, as well as the snow-dominated regime in R5.
Water 2018, 10, x FOR PEER REVIEW 12 of 22 As previously stated, these regions represent the semiarid conditions in Sierra Nevada, and their features tend to be enhanced in the future climate scenarios, especially under RCP 8.5, as expected. All resulting trends are positive (RCP 4.5: 0.006 (mm·days −1 )·year −1 and RCP 8.5: 0.01 (mm·days −1 )·year −1 ), although quite moderate under RCP 4.5. The decadal scale more clearly confirms the semiarid character of R2, with high variability, as well as the snow-dominated regime in R5.

Discussion
The analysis performed in the Sierra Nevada area provides an estimate of the long-term effects of climate evolution under two different emission scenarios in a mountain area representative of Mediterranean conditions. The results highlight how climate impacts on snow in Mediterranean regions can be significantly variable in short scale spatial domains [67][68][69] and quantify the estimated anomalies in snowfall on both an annual and decadal basis.
The two scenarios chosen for the analysis largely differ in the extent of the estimated impacts of climate on snowfall occurrence and amount on an annual and decadal basis. As expected, the more severe scenario, RCP 8.5, results in larger impacts in terms of anomalies of snowfall. This has also been reported from similar analyses in different studies across the world [14,[70][71][72]. Despite the fact that the decreasing evolution of precipitation on an annual basis shows no significant trends, the resulting snowfall regime for each scenario follows a significant decreasing trend associated with the long-term increasing trend for temperature.
This general significant increasing trend of temperature in future climate scenarios has been reported in climate impact analyses in other high mountain regions, like the Alps and the Himalayas [73][74][75][76], as well as in other semiarid mountain areas in Chile (the Andes) and in California (the Sierra Nevada) [77][78][79]. Similar values of the estimated trends for mean temperature have been found for the Alps (0.22 and 0.58 • C·decade −1 for RCPs 4.5 and 8.5, respectively), Sierra Nevada (0.18 and 0.46 • C·decade −1 for RCPs 4.5 and 8.5, respectively), and Chile (where an increase of 3 to 4 • C for the 21st century has been shown) [77,78]. The resulting trends in minimum temperature are lower, however, than those found in these other regions (e.g., 0.3 and 1.1 • C·decade −1 for RCPs 4.5 and 8.5, respectively, in the Himalayas and 0.08 and 0.26 • C·decade −1 for RCPs 4.5 and 8.5, respectively, in Sierra Nevada). This lower value is likely due to the different ranges of altitude between both sites, but scale effects associated with the downscaling in Sierra Nevada could also influence the results. In addition, maximum temperature trends show different behaviors between both areas (e.g., 0.20 and 0.30 • C·decade −1 for RCPs 4.5 and 8.5, respectively, in the Himalayas and 0.21 and 0.66 • C·decade −1 for RCPs 4.5 and 8.5, respectively, in Sierra Nevada).
The lack of significance in the trend analysis of precipitation has also been reported in these studies [80][81][82][83][84], with higher decreasing trend values observed during the 21st century under semiarid conditions (e.g., −9% and −19% for RCPs 4.5 and 8.5, respectively, in the Chilean Andes and −1.5% and −5.5% for RCPs 4.5 and 8.5, respectively, in Sierra Nevada). These results also show how apparent decreases of snowfall are expected in all future scenarios [74,76,80,81,[85][86][87], with a higher loss of snowfall amount in areas above 1500 m a.s.l. but a larger impact in lower areas were snowfall might not even occur in the future [74,76,86]. The resulting impacts on snowfall differs largely among sites. For example, decreasing trends of 32 and 65 mm·decade −1 were found in the Alps for each scenario [75], which were larger than those found in Sierra Nevada. However, in relative values, these results are closer to the decadal decrease of snowfall estimated in Canada [85], which are 5 and 10% of the current decadal values for RCPs 4.5 and 8.5, respectively. This is in line with the 2.8 and 6.1% obtained for the same scenarios in Sierra Nevada. A key result from this analysis is the likely impact of climate scenarios on the torrentiality of snowfall. Despite the mean annual trend being a decrease (with high significance) during the study period, the extreme annual values of snowfall within each decade show an extremely high variability, with large amplitudes under both scenarios in most regions. From Tables 1 and 2, it can be observed how the number of snowfall events in the future period is lower than that during the reference period, while the mean duration of each event is also lower. This behavior has also been reported in the Alps [76], with the mean value of the number of days of snowfall in future scenarios being 50% that of the reference period. That is in the range of the 33% and 76% found in Sierra Nevada under scenarios RCP 4.5 and RCP 8.5, respectively. This impact is highlighted by the resulting Isnowfall (Figure 9), which increases during the second half of the study period and reaches the highest values in regions R1 and R2, which represent the semiarid environments in Sierra Nevada. This can also be observed from the results in Table 3, which shows the proportion between the mean values for each variable over the whole future time period under scenarios RCP 8.5 and RCP 4.5 for each region in Sierra Nevada and the whole study area. Whereas the ratios of the maximum and mean temperature and precipitation are very similar among regions, the ratio of the minimum temperature shows clear differences that are associated with the frequency of snowfall in each region. This greatly impacts the snowfall under each scenario, with ratios that differ among regions for both the maximum and minimum snowfall and the associated trend, which is amplified in the more severe scenario. This is especially true for the most snow-influenced regions (R5, R4). The apparent stabilization of the snowfall anomalies during the last decades of the 21st century under RCP 4.5, compared to the continuous fall estimated under RCP 8.5, stresses the importance of mitigation actions against the adaptation actions in the climate impacts simulated by the model.  Figure 11 sums up the regional differences found in the study area for both the mean annual values of the targeted variables and their mean annual trends. It can be easily observed how temperature trends are spatially similar despite the fact that the mean original values are not homogeneous in space, whereas the precipitation pattern is variable both for the mean values and their trends. The resulting snowfall distribution is hence variable, both as mean and trend values.
The limitations of the results discussed above are mainly associated with the skill of the downscaled data and with the simplifying hypothesis in the snowfall/rainfall partition algorithm applied in this work. Additionally, new weather stations in the higher areas are required to further validate the spatial interpolation of both temperature and precipitation. Nevertheless, the results are aligned with those projected during the 21st century in other high mountain regions and provide the first distributed assessment of future climate scenario impacts on snowfall in Sierra Nevada, which has similarities with other semiarid mountain regions that allow some transferability to other less-monitored areas.

Conclusions
This work assessed climate impacts on the snowfall regime in a Mediterranean mountain area in southern Spain (i.e., the Sierra Nevada range) under two different future scenarios (RCP 4.5 and RCP 8.5). The results indicate a global decrease of annual snowfall, with a significant trend that ranges from 0.21 to 0.55 (mm·year −1 )·year −1 and changes depending on the scenario and region in the area. However, the major impact of climate on the snowfall regime is reflected in the torrential character of snowfall occurrence, with a decrease in the number of days with snowfall and a significant increase in the global mean snowfall intensity under both scenarios. This torrentiality is more extremely increased in the semiarid region of Andarax in Sierra Nevada, the area currently less influenced by snow. This may have a relevant impact on the fluvial regime, which is currently characterized by non-perennial flows, with the occurrence of flash-flood events.
Climate projections are not forecasts, but rather likely long-term pathways that climate variables may follow as they are driven by the long-term dynamics of the global circulation system in the Earth under different emission scenarios. The uncertainties associated with the estimation of snowfall occurrence from the assumptions made in this work add up to uncertainties in the global models, the downscaling techniques, the spatial interpolation algorithms, and the observational datasets. Nevertheless, the thus-projected trends of snowfall in the Sierra Nevada area show how climate impacts on the snow regime in Mediterranean and other semiarid regions are highly dependent on the trend towards torrentiality observed in the precipitation regime in these areas. These trends focus on the extreme value regimes of both precipitation and temperature variables, rather than on the mean value regime. Moreover, the spatial analysis' results highlight the large heterogeneity of these impacts on the local scale and point out the need for high resolution modeling

Conclusions
This work assessed climate impacts on the snowfall regime in a Mediterranean mountain area in southern Spain (i.e., the Sierra Nevada range) under two different future scenarios (RCP 4.5 and RCP 8.5). The results indicate a global decrease of annual snowfall, with a significant trend that ranges from 0.21 to 0.55 (mm·year −1 )·year −1 and changes depending on the scenario and region in the area. However, the major impact of climate on the snowfall regime is reflected in the torrential character of snowfall occurrence, with a decrease in the number of days with snowfall and a significant increase in the global mean snowfall intensity under both scenarios. This torrentiality is more extremely increased in the semiarid region of Andarax in Sierra Nevada, the area currently less influenced by snow. This may have a relevant impact on the fluvial regime, which is currently characterized by non-perennial flows, with the occurrence of flash-flood events.
Climate projections are not forecasts, but rather likely long-term pathways that climate variables may follow as they are driven by the long-term dynamics of the global circulation system in the Earth under different emission scenarios. The uncertainties associated with the estimation of snowfall occurrence from the assumptions made in this work add up to uncertainties in the global models, the downscaling techniques, the spatial interpolation algorithms, and the observational datasets. Nevertheless, the thus-projected trends of snowfall in the Sierra Nevada area show how climate impacts on the snow regime in Mediterranean and other semiarid regions are highly dependent on the trend towards torrentiality observed in the precipitation regime in these areas. These trends focus on the extreme value regimes of both precipitation and temperature variables, rather than on the mean value regime. Moreover, the spatial analysis' results highlight the large heterogeneity of these impacts on the local scale and point out the need for high resolution modeling to derive further information linked to the snow regime, such as runoff and river flow, soil wetness, vegetation distribution patterns, erosion rates, water temperature, and groundwater recharge.
Author Contributions: M.J.P.-P. and R.P. were responsible for processing the data. M.J.P.-P was responsible for analyzing the data and communicating with the journal. R.P. contributed technical details and analysis support and M.J.P. was responsible for overseeing the research and providing critical insight and recommendations regarding the focus, structure, and content of the paper. All authors participated in writing and proofreading throughout the publication process.