A Decrease in the Regulatory Effect of Snow-Related Phenomena in Spanish Mountain Areas Due to Climate Change

Climate change undoubtedly will affect snow events as temperature and precipitation are expected to change in the future. Spanish mountains are especially affected by that situation, since snow storage is there focussed on very specific periods of the hydrological year and plays a very important role in the management of water resources. In this study, an analysis of the behaviour of the complex snow-related phenomena in the four main mountain regions of Spain in the next 50 years is conducted. The ASTER hydrological model is applied using temperature and precipitation data as basic input, estimated under a climate change scenario. Results show different changes in the maximum and average expected flows, depending on the very different magnitude and sign of changes in precipitation. An increase of flooding episodes may occur as a result of a complex relation between changes in precipitation and an increase in maximum snowmelt intensities that range from 2.1% in the Pyrenees to 7.4% in the Cantabrian Mountains. However, common patterns are shown in a shorter duration of the snow bulk reserves, expected to occur 45 days earlier for the Cantabrian Mountains, and about 30 days for the rest of the studied mountain regions. Changes observed also lead to a concerning decrease in the regulatory effect of the snow-related phenomena in the Spanish rivers, with a decrease in the average snow accumulation that ranges from about 28% for the Pyrenees and Sierra Nevada to 42% for the Central System and the Cantabrian Mountains. A decrease in average flow is expected, fluctuating from 2.4% in the Pyrenees to 7.3% in Cantabrian Mountains, only increasing in the Central System by 4.0%, making all necessary to develop new adaptation measures to climate change.


Introduction
In the current era, where society demands sustainable development together with a respectful treatment of the environment, the availability of water resources is an essential issue both qualitatively and quantitatively. The quantification of water resources in the form of snow is therefore of high interest as it corresponds to a natural system for regulating the flow of water. Researches showed the connection between the timing and volume of snow accumulation and snowmelt in hydrological systems [1][2][3]. Snow melting creates a set of water resources which need to be controlled and quantified, establishing their flow contribution to the rivers. Forecasting both variables, snowmelt and flow contribution, is important for the ordinary management that each Basin Organization performs and acquires great relevance for managing extreme hydrological phenomena, i.e., droughts and flooding, especially when trying to envisage and mitigate the damage they may cause [4][5][6][7][8].
In Spain, the meteorological conditions of mountainous areas, as well as their geographical and geomorphological characteristics, cause the snow-related phenomena to be focussed on very specific periods of the hydrological year, mainly between December and May, while presenting great variability in time and location [9]. Accumulated precipitations in the form of snow take place in these areas between November and April, while snowmelt is centred in spring, in the March-June period. This aspect implies an important regulatory effect of water resources [10], as the incorporation of the snowfalls into the fluvial channels are delayed during winter (at that time, snow is accumulated in the mountains).
Snowmelt in the highest parts of Spanish basins starts in spring and accelerates between the end of March and April, even continuing until July in the highest areas of the Pyrenees [9,10]. Normally this melted snow is gradually incorporated into the river networks and the reservoirs as runoff (ordinary regulation). Therefore, snow regulates the water resources of the country as it ensures a base flow in the rivers during the beginning of the summer period, which is characterized by low precipitation, high evapotranspiration and significant water demand. However, in certain specific weather conditions, such as a rapid increase in temperature combined with rains, significant melting of snow in short periods may occur, leading to high flows in river channels and causing flooding [4].
In the last decades various studies have been carried out to analyze the phenomenon of snow in Spain and its contribution to flooding [10][11][12][13][14].
In this paper, the expected behaviour of the snow-related phenomena in Spain in a future scenario considering the effect of climate change and their consequence in the regulatory effect role that snow has in Spanish rivers are assessed. The main variables that control the snow-related phenomena are studied: temperature, precipitation, snowfall, snow accumulation, snowmelt and flow produced. Four Spanish mountain regions are studied: the Pyrenees, Sierra Nevada, the Central System and the Cantabrian Mountains. The results of this work contribute to improve the knowledge about the effect of climate change in water resources in the form of snow in the main mountain regions of Spain, with few previous similar studies in the same region [15].
Changes in the snow patterns in mountains and ranges are undoubtedly one of the clearest indicators that climate change is taking place on our planet since snow is extremely sensitive to changes in global meteorological conditions [16][17][18][19][20][21][22]. Snow depends on the type of precipitation, i.e., in the form of snow or liquid, and water resources associated also depend on the duration of snow in mountains and the speed of snowmelt, which leads snow to disappearing as a stored water resource, being integrated into the river network. Both aspects are related among other factors to temperature [23]. An increment in temperature, as expected in climate change scenarios, will reduce precipitation in form of snow and considerably affect snowmelt by accelerating it.
In this work, variables that control the snow-related phenomena are obtained using the hydrological model ASTER [24][25][26][27], a model applied in several countries and implemented by Spanish water authorities as an accurate model in river basins with strong snow phenomena. Results obtained enable establishing both the actual state and the future scenario considering forecasts given by climate models approved by the IPCC [28].
In the near future, a more efficient consumption of water will be required to allow a sustainable life system over time. The management of water resources in a country like Spain, where those resources are scarce, will be more complex. The study on how the snow-related hydrological processes evolve in the main mountainous areas of Spain is therefore decisive to quantify its regulatory effect and establish future management plans.

Mountain Regions Studied
The Pyrenees, Sierra Nevada, the Central System and the Cantabrian Mountains were the Spanish mountain regions selected to be studied in this work. Those areas are defined by Spanish Institutions [29] as regions where the snow-related phenomena reach a special relevance. These four regions also match with areas occupied by ice during the last ice age.
The Pyrenees extends about 490 km between Spain and France and its highest altitude corresponds to the peak of Aneto with 3404 m. The set of headwaters of the main rivers of the Spanish Pyrenees were simulated, covering an area of 10,300 km 2 with an average elevation of around 1400 m. Sierra Nevada extends about 80 km to the south of Spain and its highest altitude is the Mulhacen with 3479 m. In this region there are two delimited areas: the Guadalquivir river basin and the Mediterranean basin; the simulation covers an area of 1253 km 2 with an average elevation of around 1500 m.
The Central System extends about 600 km in the central area of the Iberian Peninsula, just north of the 40th parallel, and its highest altitude corresponds to the peak of Almanzor with 2592 m. This mountain range is in fact a succession of ranges separated by valleys or mountain passes. Gredos, Guadarrama and Ayllon are some of those ranges. In this region, the simulation covers an area of 7758 km 2 with an average elevation of around 1069 m. Unlike the two previous regions, the snow-related phenomena in the Central System have unequal importance that depends on each basin since their average levels vary between 850 m and 1500 m.
The Cantabrian Mountains extend about 480 km in the north of Spain, parallel to the Cantabrian Sea (Bay of Biscay) and their highest altitude is the Torre Cerredo with 2648 m; this range is the most western one in Europe. In this region, there are two main areas: the Duero river area and the Cantabrian Sea area. In total, the simulation covers an area of 14,099 km 2 with an average elevation of around 1169 m. Snow is especially important in the cloud fronts coming from the north since this mountain range is the first orographic obstacle found by these fronts in their movement from north to south. The fronts discharge significant precipitations in the form of snow, making the Cantabrian Mountains the largest snow-covered surface of Spain. Table 1 lists the number of basins analysed in each region as well as some statistical values regarding their size along with the actual average values of temperature and precipitation (time period 2007-2018). The location of the regions in Spain and the set of basins studied in each region are displayed in Figure 1.

Variables Analysed
The main variables that control the snow-related phenomena and their contribution in terms of water resources are temperature, precipitation, snowfall, snow accumulation, snowmelt and flow produced.
Temperature fluctuation over a year affects the precipitation in form of snow as well as melting period and speed. Temperature evolution is essential to define what will be the change in the importance of snow in each region and what will be the change in the way in which this resource becomes part of the flows in rivers (melting).
Precipitation may be in form of snow or water, with the threshold temperature between these two varying from +2 • C to −2 • C according to several atmospheric conditions. Thus, precipitation refers to the total precipitation while snowfall refers to the precipitation in form of snow. Precipitations evolution will affect the snow-related phenomena, as the more total precipitation the more snowfall, although temperature also influences the latter.

Variables Analysed
The main variables that control the snow-related phenomena and their contribution in terms of water resources are temperature, precipitation, snowfall, snow accumulation, snowmelt and flow produced.
Temperature fluctuation over a year affects the precipitation in form of snow as well as melting period and speed. Temperature evolution is essential to define what will be the change in the importance of snow in each region and what will be the change in the way in which this resource becomes part of the flows in rivers (melting).
Precipitation may be in form of snow or water, with the threshold temperature between these two varying from +2 °C to −2 °C according to several atmospheric conditions. Thus, precipitation refers to the total precipitation while snowfall refers to the precipitation in form of snow. Precipitations evolution will affect the snow-related phenomena, as the more total precipitation the more snowfall, although temperature also influences the latter.
The temperatures and the precipitations (both rainfall and snowfall) are the inputs required by the ASTER model to predict the snow accumulation, the snowmelt and the flow (model outputs).
Snow accumulation refers to the water resources accumulated in the basin in form of snow. These water resources will be incorporated into the river networks (flow) by snowmelt. Usually, snow starts to melt practically at the same time as the first snowfall begins, especially in the lower areas of the basin during the first months.
Later, with the increment of temperatures in spring, melting starts in the highest areas where snow was accumulated throughout the winter. This leads to a rapid increase in the contributions (flow) due to this phenomenon. Thus, the alteration in the temperature The temperatures and the precipitations (both rainfall and snowfall) are the inputs required by the ASTER model to predict the snow accumulation, the snowmelt and the flow (model outputs).
Snow accumulation refers to the water resources accumulated in the basin in form of snow. These water resources will be incorporated into the river networks (flow) by snowmelt. Usually, snow starts to melt practically at the same time as the first snowfall begins, especially in the lower areas of the basin during the first months.
Later, with the increment of temperatures in spring, melting starts in the highest areas where snow was accumulated throughout the winter. This leads to a rapid increase in the contributions (flow) due to this phenomenon. Thus, the alteration in the temperature regime could cause melting in a basin to be advanced, making it coincide with raining periods and increasing the flows of liquid water in the rivers, raising the risk of flooding and affecting the current management of a basin.

ASTER Hydrological Model
ASTER is a continuous simulation distributed hydrological model that simulates the hydrological behaviour of a high mountain sub-basin and calculates the expected stream flow in a river and the equivalent volume of water stored at any time in the studied basin using precipitation and temperature as the main inputs variables. This model, used currently by Spanish water authorities [25,26] was specially developed for basins with steep reliefs, with very marked meteorological changes in time and space [30,31]. Some remarkable features of the ASTER model include the possibility of choosing different spatial and temporal series for the simulation, the fragmented treatment of precipitation and temperature variables and the specific and distributed analysis of the snow-related phenomena. In this work a daily time resolution was used, matching the time series resolution provided by the climate models. The daily accumulated precipitation and the average daily temperature were the input variables of the model, from which the evolution of the snow-related phenomena and the daily flows were obtained.
ASTER model calculation algorithm follows the classical hydrological models based on quasi-Distributed continuous simulation deposits similar to the Canadian CEQUEAU model [32], with a variable temporal operating scale. The snowmelt routine is based on the studies carried out by Anderson [33] and applied in the US model NWSRFS of the National Weather Service [34][35][36]. Particularly, the model calculates snowmelt and snow accumulation making use of an energy balance equation [25,26].
For defining the snow accumulation, quantification of snowfall is conducted based on a given temperature named "rain/snow temperature" T rain/snow which is defined for each basin, a typical value being around 1.5 • C. For that value, half precipitation is considered to fall as snow and half as water. If precipitation temperature is greater than T rain/snow by 2 • C, all the precipitation is considered as liquid form. If precipitation temperature is lower than T rain/snow by 2 • C, all the precipitation is considered to fall as snow. A linear interpolation is made between the previous values.
Snowmelt is equal to the result of two components: snowmelt due to precipitation (M p ) and snowmelt without rainfall (M c ).
The snowmelt due to precipitation (M p ) is the snowmelt produced by the precipitation energy and it is calculated using Equation (1), assuming the snow surface temperature is equal to 0 • C [33]: where P is the precipitation, in mm; T p is the temperature of the precipitation, in • C; and f p is a coefficient that estimates the fraction of precipitation in the form of rain, its value being set to 0.5 for T rain/snow . Besides the precipitation energy, other secondary terms like radiation and condensation energy may be added to calculate snowmelt when precipitation occurs. The snowmelt without rainfall (M c ) corresponds to the snowmelt on a dry day and is calculated using Equation (2), which is based on non-forced convection (note that this kind of estimations are commonly used in hydrological models along with other empirical relations depending on air temperature since estimating wind speed is difficult) [35]: where T a is the air temperature, in • C; T mi is the snowmelt temperature in the basin, in • C; and M f is the melt factor, in mm/ • C, which varies seasonally and depends on several implicit factors, mainly on radiation (latitude, orientation, slope, albedo and vegetation cover), exposition and wind. This factor may be expressed by Equation (3) where n is the number of days from 21 March. It should be noted that the value T mi (snowmelt temperature) may be different from 0 • C. For instance, on calm and clear days, when solar radiation dominates the energy balance, melting can occur below 0 • C, while on clear nights, when the long wave radiation emitted is significant, there would be no melting until air temperature exceeds 0 • C. Therefore, T mi needs to be calibrated in each case analysed. By using the previous equations, the amount of water incorporated in the river networks can be obtained and added to the one directly produced by rain (precipitations in form of water) to estimate the flow expected on each period of the year.

Selection of the Climate Model
The input data required by the ASTER hydrological model are the temperatures and the precipitations (see Section 2.2), stemming from climate models.
In this work, temperature and precipitations time series provided by the Spanish Agency of Meteorology AEMET [29] are considered. In particular, among the 24 climate models made available by AEMET, ANA-MRI-CGCM3-rcp85 was selected for conducting the present study. This model provides regionalized climate projections from the CMIP5 (Coupled Model Intercomparison Project phase 5), according to the IPCC and framed within the greenhouse gas emission scenario (RCP 8.5), by using an analogic downscaling method (ANA).
The motivation for using the ANA-MRI-CGCM3-rcp85 model results is twofold: (i) previous experience [4]; and (ii) current validation ( Figure 2). Previous experience refers to a previous work carried out by the authors [4], where the analyses showed that this climate model achieved a good fit with real annual-average daily-flow-rates. Besides, before using the model again for this study, the authors checked and demonstrated, for the particular case of the Arga river basin (see Figure 2), that the validity of the ANA-MRI-CGCM3-rcp85 model predictions still stands.
The input data required by the ASTER hydrological model are the temperatures a the precipitations (see Section 2.2), stemming from climate models.
In this work, temperature and precipitations time series provided by the Span Agency of Meteorology AEMET [30] are considered. In particular, among the 24 clim models made available by AEMET, ANA-MRI-CGCM3-rcp85 was selected for conducti the present study. This model provides regionalized climate projections from the CMI (Coupled Model Intercomparison Project phase 5), according to the IPCC and fram within the greenhouse gas emission scenario (RCP 8.5), by using an analogic downscali method (ANA).
The motivation for using the ANA-MRI-CGCM3-rcp85 model results is twofold: previous experience [4]; and (ii) current validation ( Figure 2). Previous experience ref to a previous work carried out by the authors [4], where the analyses showed that t climate model achieved a good fit with real annual-average daily-flow-rates. Besides, fore using the model again for this study, the authors checked and demonstrated, for t particular case of the Arga river basin (see Figure 2), that the validity of the ANA-M CGCM3-rcp85 model predictions still stands. Figure 2 shows the results of a thorough simulation study on the extreme rainf predictions of several climate models (provided by AEMET), over the period 2011-20 on the Arga river basin. In the figure, the results of the numerical study are compar with experimental data measured by the local authorities. According to Ebro Water A thority, Arga river basin is one of the best monitored basins in the region, being the perimental data reliable and thus suitable for validation purposes. From Figure 2, it c be concluded that ANA-MRI-CGCM3-rcp85 is the most accurate model and it is thus us as reference climate model for this work. Precipitación (mm/año) Figure 2. Estimated precipitation in the Arga river basin using several climate models provided by AEMET [28] for the time period 2011-2017 against the real observations measured by local authorities at that period. Figure 2 shows the results of a thorough simulation study on the extreme rainfall predictions of several climate models (provided by AEMET), over the period 2011-2017 on the Arga river basin. In the figure, the results of the numerical study are compared with experimental data measured by the local authorities. According to Ebro Water Authority, Arga river basin is one of the best monitored basins in the region, being the experimental data reliable and thus suitable for validation purposes. From Figure 2, it can be concluded that ANA-MRI-CGCM3-rcp85 is the most accurate model and it is thus used as reference climate model for this work.

Calibration and Validation of the ASTER Model
For each hydrological basin studied in this work, the calibration of the ASTER hydrological model is done against the flow data measured at each of the rivers defining each of the hydrological basins studied, obtained from the recorded level-flow data in gauging stations (flow control points) for the reference time period, which in this case it extended from 1 October 2006 to 30 September 2010. The model is calibrated defining the hydrological parameters that make it possible to optimize the correlation coefficient and the Nash-Sutcliffe parameter (NSE) by statistically analysing the observed flow series and those estimated by the ASTER model. The resulting model is then validated by comparing the flow model predictions from 1 October 2010 to 30 September 2018 with the corresponding available measured data.
Not only the flow, but also the snow accumulation predicted by the ASTER hydrological model in the four regions under analysis is confirmed with measured snow depth and snow density gathered from campaigns carried out three times per year at the beacon network during the last 30 years (EHRIN program) [29]. Table 2 shows a summary of the model parameters resulting from the calibration process. The calibration and validation carried out in a previous work for a longer period  in all the basins of the Pyrenees [9], showed a good consistency of the model to the variability that might be associated with climate change. However, several uncertainties still remain associated with changing climate conditions, such as infiltration due to evolution of vegetation, which may not be easy to predict and may be influenced as well by other anthropic effects non-related with climate change. As snowmelt is calculated based on temperature, other uncertainties associated with changes in incoming radiation or other parameters will result in temperature changes that will be taken into account by the model. As a distributed model, spatial interpolation and temperature and precipitation elevation gradients enable to transform grid data from climate models over complex terrain.  Table 3 shows the results of the validation of the ASTER model from 1 October 2010 to 30 September 2018, for several basins, representing the four studied mountain regions and for two different periods (the complete period and the snowmelt period). The average observed flow and the average flow that is calculated by the ASTER model are presented along with the coefficient of determination (R 2 ) and the Nash-Sutcliffe parameter (NSE), suggesting good predictive skills for the ASTER hydrological model.

Results of the Projections for the Selected Climate Model
Using the values reported in Table 2 and the temperature and precipitations given by the ANA-MRI-CGCM3-rcp85 climate model, the ASTER hydrological model was used for conducting a future projection of the snow-related phenomena in the period 2041-2070. Table 4 shows the resulting values of the temperature, precipitation, snowfall, snow accumulation, snowmelt and flow in the four regions analysed. For each parameter, the maximum and average values are listed, both in terms of the present-day scenario (from 2006 to 2018) and the future scenario, along with its relative change (Rel. change) as a percentage following Equation (4): where V act is the actual value of a parameter and V f ut is the future estimated value.  Maximum values represent the highest daily values found in the present-day scenario (from 2006 to 2018). Those maximum values cannot be directly compared with the ones corresponding to the future scenarios: the difference in the length of both periods may lead the maxima in the future scenario to representing much more rare events (once in 30 years) than those during the present-day period (once in 13 years). Consequently, a homogenization procedure was conducted for the future scenario results consisting of considering only the last 13 years of the climate projection (from 2058 to 2070) for computing the maximum values (highest daily values found in that period). Average values were not affected by that procedure and were computed for the whole 30-year period (from 2041 to 2070).

Pyrenees
Temperatures in the Pyrenees show a great increase in annual average values, with a moderate increase in the maximum value. Precipitations are expected to show a slight increase in terms of average values but a decrease in the maximum values. Both average and maximum snowfall values decrease due to the increment in average temperatures and the decrease of maximum precipitations. Figure 3 shows average values for precipitation, accumulated snow, snowmelt and flow in the Pyrenees. Figure 3a confirms the slight increase in precipitations and it also reveals a change in the precipitation pattern, decreasing in summer-autumn and increasing in winter and spring. Figure 3b shows that the average accumulated snow decreases throughout the year. Even though the average snowmelt decreases for the whole year, the increase in temperatures results in an increase in snow meltdowns from January to April as displayed in Figure 3c. This also means a change in the snowmelt pattern, which will be forwarded in time and will take place more than a month before to the present-day scenario. Average flow decreases along the year too, but the increase in precipitations and snowmelt during the mentioned period leads to an increase in the flow stream to rivers during the first months of the year (Figure 3d). Advance in snowmelt and the decrease in precipitations in summer also result in decreasing the average flow between the months of June and October. as displayed in Figure 3c. This also means a change in the snowmelt pattern, which will be forwarded in time and will take place more than a month before to the present-day scenario. Average flow decreases along the year too, but the increase in precipitations and snowmelt during the mentioned period leads to an increase in the flow stream to rivers during the first months of the year (Figure 3d). Advance in snowmelt and the decrease in precipitations in summer also result in decreasing the average flow between the months of June and October.

Sierra Nevada
Temperatures in Sierra Nevada show a great increase in annual average values, while a slight decrease is observed in the maximum values. Precipitations increase, the increment being moderate for both the average and the maximum values. Maximum snowfall slightly increases in line with the total precipitation, although the average snowfall decreases due to the increment in average temperatures. Figure 5 shows average values for precipitation, accumulated snow, snowmelt and flow in Sierra Nevada. Figure 5a displays how precipitations are concentrated in autumn in the present-day scenario, while in the future scenario precipitations will occur mainly in winter and early spring. This will cause a decrease in the accumulated precipitations in summer and autumn and reveals a change in the precipitation pattern. Average accumulated snow (Figure 5b) decreases throughout the year. Snowmelt average value decreases for the whole year, but Figure 5c shows a clear increase in average snowmelt from January to April, also leading to an advance in the snowmelt pattern (in the actual time, snowmelt mainly occurs in spring) due to the increase in temperatures similar to what was observed in the Pyrenees. A decrease in average annual flow is expected but also showing a pattern change: flow will decrease from June to January and increase in late winter and spring ( Figure 5d).

Central System
Temperatures in the Central System show a great increase in average values and a slight decrease in maximum values. Similarly, a moderate increase in average precipitations is expected while maximum values will show a decrease. Snowfall decreases both in terms of maximum and average values due to the increment in average temperatures. Figure 7 shows average values for precipitation, accumulated snow, snowmelt and flow in the Central System. Figure 7a reveals a change in the precipitation pattern with precipitations decreasing in autumn and increasing along the months of January to April. Consequently, snow accumulations and flow rates will experience significant changes in both magnitude and temporal location with respect to the actual values. Figure 7b shows that the average accumulated snow decreases throughout the year, while Figure 7c displays an increase in snowmelt from January to March due to the increase in temperatures. Average annual flow experiences a slight increase, especially in months between January and April caused by the increase in snowmelt (Figure 7d). Figure 8 shows the evolution of the previous variables in terms of maximum values. According to Figure 8a, maximum precipitation values in the future scenario will be higher in winter. Figure 8b shows that maximum snow accumulation will decrease throughout the whole year. Maximum snow accumulations will reach their peak in winter and will start decreasing from March, evidencing a change in the melting pattern which is also confirmed in Figure 8c. In the present-day scenario, the maximum snowmelt occurs in spring while in the future scenario is expected to take place in winter and decrease in the months of April and May. Therefore, melting will be forwarded and accelerated with respect to the actual situation. A moderate increase in the maximum flow value is detected and as Figure 8d shows in the Central System winter will favour the increase of flow due to the early melting expected to occur in those months. The advancement in the meltdown pattern together with the increase in total precipitations in the first half of the year will result in increasing the maximum flow in this period. The maximum flow is expected to decrease in autumn as a consequence of the decrease in precipitations and the increase in temperatures.
in spring while in the future scenario is expected to take place in winter and decrease in the months of April and May. Therefore, melting will be forwarded and accelerated with respect to the actual situation. A moderate increase in the maximum flow value is detected and as Figure 8d shows in the Central System winter will favour the increase of flow due to the early melting expected to occur in those months. The advancement in the meltdown pattern together with the increase in total precipitations in the first half of the year will result in increasing the maximum flow in this period. The maximum flow is expected to decrease in autumn as a consequence of the decrease in precipitations and the increase in temperatures.

Cantabrian Mountains
Temperatures in the Cantabrian Mountains show an increase in both annual average and maximum values. Precipitations show a slight decrease in average values and a great decrease in maximum values. Average snowfall decreases due to the increment in average temperatures and the decrease in average precipitations, but maximum snowfall increases, which indicates that important heavy snowfalls are expected to occur in the future scenario. Figure 9 shows average values for precipitation, accumulated snow, snowmelt and flow in the Central System. Figure 9a displays a change in the precipitation pattern similar to the ones observed in the previous regions studied, with precipitations decreasing in autumn and increasing in winter. Figure 9b confirms that the increase in temperatures and the decrease in average snowfall results in a significant decrease in the accumulated snow throughout the year. A change in the melting pattern is also observed in Figure 9c, as the highest average snowmelt values move from spring (present-day situation) to winter (future situation). Therefore, melting will be forwarded and accelerated with respect to the actual situation. A shorter duration of the snow mantle is also evidenced (difference length of the "plateaus" observed in Figure 9c). Figure 9d Figure 10 shows the evolution of the previous variables in terms of maximum values. Figure 10a confirms the decrease in maximum precipitation. Figure 10b shows that the maximum accumulated snow also decreases throughout the whole year. The large increase in the maximum snowfall but the decrease in accumulated snow, points to a change in the snow pattern in which heavy snowfall will be favoured but their permanence on the peaks will be shorter due to the significant increase in temperatures. According to Figure 10c, the maximum snowmelt is expected to increase between January and April and decrease along the months of April and May. That confirms the change in the melting pattern observed with average values. As both the maximum snowmelt and maximum precipitation decrease, the maximum flow also decreases. Figure 10d shows how the flow will increase during the first months of the year as a result of snowmelt, despite the slight decrease in precipitations, which evidences that the maximum flow will be influenced by snowmelt in a greater proportion than in the present time.  Figure 10 shows the evolution of the previous variables in terms of maximum values. Figure 10a confirms the decrease in maximum precipitation. Figure 10b shows that the maximum accumulated snow also decreases throughout the whole year. The large increase in the maximum snowfall but the decrease in accumulated snow, points to a change in the snow pattern in which heavy snowfall will be favoured but their permanence on the peaks will be shorter due to the significant increase in temperatures. According to Figure 10c, the maximum snowmelt is expected to increase between January and April and decrease along the months of April and May. That confirms the change in the melting pattern observed with average values. As both the maximum snowmelt and maximum precipitation decrease, the maximum flow also decreases. Figure 10d shows how the flow will increase during the first months of the year as a result of snowmelt, despite the slight decrease in precipitations, which evidences that the maximum flow will be influenced by snowmelt in a greater proportion than in the present time.

Radiation Influence
ASTER hydrological model has been calibrated and validated in the Pyrenees for a long period that ranges from 1950 to 2005 [9], showing a good fitting for the already changing climate variables caused by climate change that were non-explicit, such as radiation. However, greater changes in radiation are expected in the next 50 years, so the model might show uncertainties associated with this variable. Radiation depends on multiple factors that are difficult to predict. Besides, other effects such as snow impurities affecting snow albedo may have a greater effect on melting than the radiation input itself [38].
Nevertheless, melting related with radiation is not very significant in Spanish mountains, as both average and maximum snowmelt take place during medium radiation months. This is different to what occurs in other latitudes or places with higher elevations, where snow water equivalent storage is still important during high radiation periods.
Other limitations of the model that may lead to uncertainties have already been stated in a previous section, such as infiltration due to future development of vegetation, which may not be easy to predict and may be influenced as well by other anthropic effects.

Decrease in the Regulatory Effect and Increase of Extreme Events
The previous analyses evidence a tendency to a change in the pattern of the snowrelated phenomena in Spain, which will decrease its regulatory effect on rivers. Nowadays the influence of snow on the contribution to water resources in Spain is highly variable in each of the different regions analysed. For instance, and according to the work of Arenillas et al. [9] who studied snow-related phenomena in Spain during the period from 1950 to 2007, in the Central System, the snow at its highest point (spring) represents a total of 10% of the rainfall received but this value reduces to 5% for the entire hydrological year. In regions where snow behaviour is much more accentuated, such as Sierra Nevada or the Pyrenees, these values can rise to a maximum of 55% and 40% respectively, leaving values

Radiation Influence
ASTER hydrological model has been calibrated and validated in the Pyrenees for a long period that ranges from 1950 to 2005 [9], showing a good fitting for the already changing climate variables caused by climate change that were non-explicit, such as radiation. However, greater changes in radiation are expected in the next 50 years, so the model might show uncertainties associated with this variable. Radiation depends on multiple factors that are difficult to predict. Besides, other effects such as snow impurities affecting snow albedo may have a greater effect on melting than the radiation input itself [37].
Nevertheless, melting related with radiation is not very significant in Spanish mountains, as both average and maximum snowmelt take place during medium radiation months. This is different to what occurs in other latitudes or places with higher elevations, where snow water equivalent storage is still important during high radiation periods.
Other limitations of the model that may lead to uncertainties have already been stated in a previous section, such as infiltration due to future development of vegetation, which may not be easy to predict and may be influenced as well by other anthropic effects.

Decrease in the Regulatory Effect and Increase of Extreme Events
The previous analyses evidence a tendency to a change in the pattern of the snowrelated phenomena in Spain, which will decrease its regulatory effect on rivers. Nowadays the influence of snow on the contribution to water resources in Spain is highly variable in each of the different regions analysed. For instance, and according to the work of Arenillas et al. [9] who studied snow-related phenomena in Spain during the period from 1950 to 2007, in the Central System, the snow at its highest point (spring) represents a total of 10% of the rainfall received but this value reduces to 5% for the entire hydrological year. In regions where snow behaviour is much more accentuated, such as Sierra Nevada or the Pyrenees, these values can rise to a maximum of 55% and 40% respectively, leaving values of 45% and 28% for the entire hydrological year respectively. It should be noted that higher elevation, as is the case of Sierra Nevada and the Pyrenees, implies longer snow permanence and greater regulatory effect.
Melted snow streams to rivers, increasing their flow significantly. Therefore, the change in melting patterns also leads to a change in flow patterns. As melting will move from spring to winter, maximum flows are also expected to significantly increase during winter. The hypothesis of slower snowmelt in a warmer world due to contraction of the snowmelt season to a time of lower available energy [38], must be combined with the response of the snowpack to the magnitude and sign of changes in precipitation [39].
In the Pyrenees, the water accumulated in the form of snow in the mountains for the climate change scenario undergoes an early meltdown of about 30 days compared to the current situation. This implies a lower regulatory effect from the point of view of water resources. The flow increases in the rivers during the winter season and at the beginning of spring, while the contribution of flows from the month of April is reduced. This scenario denotes a lower regulatory effect of circulating flows and therefore will force different management of the reservoirs to the current one, requiring greater regulation capacity to achieve supply conditions similar to the actual ones. A decrease on snow accumulation about 20% for every 1 • C, and a reduction on the duration of the snowpack about 20-30 days per • C was estimated for a small snow-dominated basin in the Spanish Pyrenees [39].
In Sierra Nevada melting of the accumulated snow is advanced in around 30 days. This advance is directly reflected in flowing flows, which suffer a decrease in the months of greatest demand. Therefore, the future scenario shows that the regulatory effect of water resources associated with the snow-related phenomena will be reduced. In the Central System, climate change will lead to snowmelt increasing until reaching its maximum value in February, approximately 30 days earlier compared with the time it normally occurs today. This circumstance implies an increase in circulating flows between the months of January and March. In the Cantabrian Mountains, a reduction of the regulatory effect of the accumulated snow in mountain areas is also detected, producing an earlier snowmelt of about 45 days and a much more accentuated low water level in rivers from April. These results agree with previous studies predicting that by the end of the twenty-first century, the peak flow of many North American and European rivers is expected to occur 30-40 days earlier as a result of earlier snowmelt [40,41], but changes in precipitation caused by climate change can increase or reduce flow rates.
All in all, the effect of climate change on the snow-related phenomena shows that the regulatory effect of the snow with respect to the incorporation of runoff into the riverbeds will be less significant in the future. This will require different management and probably force a greater need for reservoir capacity to achieve optimal regulation and maintain current supply guarantee conditions for the populations and agricultural lands. At the same time, extreme flood situations will continue to occur, requiring new flood risk management plans.
Especially, in terms of future flood events, the change in the snow behaviour pattern gives rise to increasing the risk of flooding and their damaging effects, with an expected increase in the number of episodes with strong flow peaks. Future snowmelt time in the four regions analysed will coincide with the first half of the year, where precipitations are also expected to increase. Besides, an increase in the average precipitations is expected in the Pyrenees, Sierra Nevada and the Central System. In the Cantabrian Mountains, precipitations along the year may decrease, but the heavy snowfall episodes expected along with the shorter duration of snow on the peaks will also cause an increased risk of flooding.
It is important to mention that all the previous analyses are entirely based on a single climate change scenario, so even though the climate model used (ANA-MRI-CGCM3-rcp85) showed a good match with the present state in Spain, some uncertainty may arise. Previous works have shown a decrease in the magnitude of extreme floods for climate model projections in Spain [36]. This is partly due to the use of global climate models from the Spanish Meteorological Agency (AEMET), which tend to provide lower extremes with a smaller variability [42]. However, it must be considered that snow accumulated is very sensitive to weather conditions, and a generalized increase in temperatures will cause a rapid and important melting and a higher probability of direct rainfall instead of snow accumulation that will provide significant water resources added to those already streaming. Thus, climate change scenarios show a greater risk of flooding situations with more damaging effects. Particularly, the increase in temperatures and decrease in average accumulations, or the delay in the accumulation of snow, may generate extreme events in which all the precipitation occur in liquid form (water), instead of part of the same being retained in the form of snow in the highest parts of mountains. This phenomenon has already happened in December 2019 in the town of Reinosa [4], which suffered the worst flooding in memory (return period estimated about 300 years) causing great damage and affecting residents.

Conclusions
This paper has studied the expected behaviour of the complex snow-related phenomena in Spain in the next 50 years according to a climate change scenario. The results obtained show a great variability for the studied mountain regions and have made it possible to quantify the variation in the behaviour of snow resources and how this will affect the circulating flows through the rivers. The main variables that control the snow-related phenomena were analysed: temperature, precipitation, snowfall, snow accumulation, snowmelt and flow produced. The most important aspects related to the snow-related phenomena were assessed, such as changes in the periods of snow, in the expected maximum and average values, duration of the bulk of the snow reserves and intensity and duration of the snowmelt process. All these changes are essential to establish their influence and contribution in flooding as well as for the optimal management of water resources, considering the regulatory effect that snow-related phenomena have in Spain (its effective contribution is delayed in time).
From a flooding risk point of view, climate change will lead to more flooding scenarios due to both a decrease of the snowfall that ranges in average from 6.5% to 15.7% and is not related to a decrease in the average precipitation, and very rapid meltdowns associated with increases in temperature that may be accompanied by liquid rainfall, revealing a future greater frequency of flooding and with higher intensity than nowadays. However, these changes in the maximum flows result from a complex relation between changes in the snow-related phenomena and changes in precipitation. As a matter of fact, an increase in maximum snowmelt intensities of 2.1% in the Pyrenees and 7.4% in the Cantabrian Mountains, lead to a respective increase of 8.0 % and 2.0% for maximum flows, while a decrease in maximum snowmelt intensities that range from 15.7% in the Cantabrian Mountains to 12.0% in Sierra Nevada, lead to a maximum flow decrease of 7.1% for the first, but an increase of 58.9% for the latter.
More concerning are the implications of the climate change scenario in the regulatory effect of the water resources, with a decrease in average flow fluctuating from 2.4% in the Pyrenees to 7.3% in Cantabrian Mountains, only increasing in the Central System by 4.0%. Therefore, this work predicts results that may be very important to implement future water policies. Results show that in a future scenario the volume of water stored in the form of snow in the upper areas of the basins will generally be significantly reduced, about 42% in Central System and Sierra Nevada and about 28% in the Pyrenees and Sierra Nevada, while average snow melting will be initiated approximately 45 days earlier than presently for the Cantabrian Mountains and 30 days for the rest of the studied mountain regions. All these issues will cause premature incorporation of runoff from snow melting, reducing up the average snowmelt from 6.6% in the Pyrenees to 15.7% in Cantabrian Mountains, and leading to a greater flow circulating in winter months and a reduction of flows in spring and summer.
However, the greater demand for water resources will continue to occur in spring and summer, which will create a greater mismatch in terms of the availability of water resources and the demand for them. This future scenario is very worrying since it will require two essential modes of action to maintain the guarantee of supply based on sustainability criteria. In the first place, a demand adjustment will be required, i.e., a more sustainable use of water resources, consumption adjustment and reuse. Secondly, more efficient management of the existing regulation reservoirs will be needed in accordance with the new climate change scenario. If these two lines of action are not developed, a serious gap will be created between demand and available resources, jeopardizing the current balance.
However, the previous analyses are entirely based on the results of a simplified distributed hydrological model, on a single climate change scenario and for a single climate model, thus some uncertainty may arise. Further investigation is therefore needed to confirm the conclusions of this work.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy reasons.