Extreme values of snow-related variables in Mediterranean regions: trends and long-term forecasting in Sierra Nevada (Spain)

Mountain areas in Mediterranean regions constitute key monitoring points for climate variability and its impacts, but long time datasets are not always available due to the difficult access to high areas, relevant for capturing temperature and precipitation regimes, and the predominance of cloudy remote sensing images during the snow season. Sierra Nevada National Park (South Spain), with altitudes higher than 3500 m a.s.l., is part of the Global Change in Mountain Regions network. Snow occurrence just 40 km from the seaside determines a wide range of biodiversity, a snowmelt fluvial regime, and the associated ecosystem services. This work presents the local trend analysis of weather variables at this area together with additional snow-related variables. For this, long term point and distributed observations from weather stations and remote sensing sources were studied and used as input and calibration datasets of a physically based snow model to derive long term series of mean and maximum daily fraction of snow covered area, annual number of days with snow, annual number of days with precipitation, mean and maximum mean daily snow water equivalent, and snowmelt and evaporation volumes. The joint analysis of weather and snow variables showed a decrease trend in the persistence and extent of the snow cover area. The precipitation regime, rather than the temperature trend, seems to be the most relevant driver on the snow regime forcing in Mediterranean areas. This poses a constraint for rigorous scenario analysis in these regions, since the precipitation pattern is poorly approximated by climatic models in these regions.


Introduction
Climate variability impact on the hydrological regime is more evident over mountain regions due their particular extreme conditions (Beniston, 2003;IPCC, 2007).Hence, a continuous monitoring of snow state variables can help to assess climate variability and early detection of shifts at significant time scales.For example, the decrease of snow cover area is an indicator of a warming climate phase conditions.Moreover, in regions where the typical alpine-mountain climate is modified by additional meteorological drivers(i.e. the recurrence of drought periods and torrential rainfall events), these impacts are enhanced, which makes them crucial areas to monitor these climate variations.This is the case of Mediterranean mountainous areas, where alpine and semiarid conditions coexist (Giorgi, 2006).
However, snow monitoring constitutes arduous work.On one hand, the limited accessibility to these sites during the snow season, and the hard working conditions for the instrumentation, make in situ continuous monitoring systems difficult to maintain (De Walle and Rango, 2008).On the other hand, the potential use of satellite remote sensing information (i.e.Landsat TM, ETM+ and OLI with temporal resolution of 16 days and spatial resolution of 30 m × 30 m, MODIS, daily images with 250 m × 250 m spatial resolution, and NOAA, 1 km×1 km daily images) to capture snow evolution at medium and large scales is limited to recent decades, poses a constraint for snow applications in heterogeneous ar-Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.
eas due to their fixed and not always high enough spatial resolution, and is not feasible during a significant fraction of the snow season because of the recurrent presence of cloud cover.Therefore, a lack of detailed snow information is common over mountainous areas, especially at high altitudes.
The use of physically based snowmelt-accumulation model, locally calibrated and validated by means of remotely sensed data, allows estimating the snowpack evolution and its significant state variables during these gap periods, over the non-monitored areas, and beyond the snow monitoring period (Pimentel et al., 2015).For that, a correct representation of the meteorological forcings, the main inputs in snow modelling, is crucial to obtain accurate snow simulations, the weather data series available in the study area, and their quality, being a relevant limitation for a satisfactory model performance.Weather variables long-term series can, thus, be used to generate long-term snow-related variables series, which allows the understanding of the past snow dynamics, and its uncertainty, the estimation of future trends, and the forecasting of snow behaviour under given climatic scenarios.The objective of this work is to study the meteorological patterns in Sierra Nevada Mountain, an alpine climate region in Southern Spain, in relation to selected snow-related variables to obtain long-term trends in this area as a basis to assess future climate scenarios.Special emphasis has been done in the analysis of the occurrence of extreme events and its impact on snow dynamics.

Study site and available data
Sierra Nevada Mountains, in Southern Spain, are a linear mountainous range parallel to the shoreline of the Mediterranean Sea (Fig. 1), where the highest summits of the Iberian Peninsula can be found (Mulhacen peak, 3479 m).It is the second highest range in Europe, only surpassed by the Alps.Consequently, snow usually appears throughout the year in the areas above 2000 m.Its proximity to the sea, only 40 km south, generates a very specific climate as the result of the interaction between sea and mountain conditions.The typical alpine climate is then modified by the Mediterranean surrounding features and, thus, snow is significantly affected.Annual precipitation fluctuates widely and can range from 300 to 1000 mm, with a high spatial variability throughout the area due to topographic effects.The average temperature ranges from −5 to 5 • C during the snow season, although minimum values of −20 • C can be found at certain times in winter.Due to these particular conditions, it belongs to the Global Change in Mountain Regions network (GLOCHAMORE) and it is recognizes as both Natural and National Parks.
Different weather station networks are available in this area, lacking, in general, data over 2000 m a.s.l.; in this work, validated daily datasets of precipitation, temperature, wind velocity, and relative humidity from selected stations of the Spanish Meteorological Agency (AEMET) were employed as input of a distributed physically based snow model (Fig. 1), calibrated and validated in this region in previous works (Herrero et al., 2010;Pimentel et al., 2014;Pérez-Palazón et al., 2014).The reference period 1960-2000 was analysed for long-term trends and basis of future scenarios assessment.The study is limited to area in which five different regions were identified to assess the local impact of weather trends on snow (Fig. 1).

Methods
The average annual precipitation and mean, maximum and minimum average daily temperature for the reference period, 1960-2000, were obtained over the study area in Sierra Nevada, together with selected annual snow related-variables that were simulated during this period using the physical and distributed hydrological model, WiMMed (Herrero et al., 2011; http://www.ugr.es/~herrero/wimmed/),developed and calibrated in this region (Herrero et al., 2009(Herrero et al., , 2011;;Millares et al., 2009;Aguilar et al., 2010;Aguilar and Polo, 2011;Pimentel et al., 2014).The model includes specific spatial interpolation algorithms for the weather variables to include the abrupt topography of this area (Herrero et al., 2010), and an energy and water balance snow model.Long-term trends of both the meteorological and snowrelated variables were analysed emphasizing the importance of extreme precipitation events.

Snow modelling
The snowmelt-accumulation model for Mediterranean sites developed by Herrero et al. (2009) is a physical model based on a point mass and energy balance, which is extended to a distributed way by means of depletion curves (Herrero et al., 2011;Pimentel et al., 2014).
The model assumes a horizontally uniform snow cover distributed in one vertical layer.In the control volume defined by the snow column per unit area, which has the atmosphere as an upper boundary and the ground as a lower one, the water mass and energy balance can be expressed by: where SWE (snow water equivalent) is the water mass in the snow column, and u is the internal energy per unit of mass (U for total internal energy).In the mass balance, R defines the precipitation rate; E is water vapour diffusion rate (evaporation/condensation); W represents the mass transport rate due to wind; and M is the melting water rate.On the other hand, regarding energy fluxes, K is the solar or short wave radiation; L the thermal or long-wave radiation; H the exchange of sensible heat with the atmosphere rate; G the heat exchange with the soil rate; and u R , u E , u W and u M are the advective heat rate terms associated with each one of the mass fluxes involved in Eq. ( 1).This approach permits an easy extension to a distributed model by means of making calculations simultaneously in each cell.However, a direct extension cannot be done when the cell area is not completely covered by snow.In such cases, the parameterization of subgrid processes is made by including depletion curves (Luce and Tarboton, 2004;Herrero et al., 2010;Pimentel et al., 2014).
The point model was calibrated and validated by means of local in situ measurements (Herrero et al., 2009) and the distibuted model thanks to remote sensing information (Herrero et al., 2010;Pimentel et al., 2014).

Long term data analysis
Annual and decadal regimes for the reference 40-year period  were analysed for the average annual precipitation (P ) and mean, maximum and minimum average daily temperature (T mean , T max , T min ) in the area, obtained from daily measurements at the selected weather stations, and the following snow related variables obtained from the validated snow model: (1) mean and maximum daily fraction of snow covered area (SC mean and SC max ), obtained as mean and maximum value of daily percentage of snow covered area over the whole area (m 2 m −2 ); (2) annual number of days with snow (n day with snow ), calculated as the number of day in which the snow exceeds 10th percentile of SC mean distribution (days); (3) annual number of days with precipitation (n day with precipitation ), counted as days with precipitation greater than 1 mm day −1 (days); and (4) mean and maximum mean daily snow water equivalent (SWE mean and SWE max ), calculated as mean and maximum values of daily average SWE (mm).Their individual temporal trend and main statistic descriptors were performed.The significance of these temporal trends was tested by means of Mann-Kendall test (Gibbons and Chakraborti, 2010).The test assumes as null hypothesis: non-monotonic trend (H 0 ), which was tested versus the alternative hypothesis: presence of monotonic trend (H 1 ).A level of significance (α) equal to 0.05 is used to define the rejection region.The snow regime dependence on weather variables was also assessed.

Extreme data analysis
The extreme data analysis was performed on a precipitation event definition basis.A precipitation event is defined as a period in which rainfall over 1 mm day −1 is registered at some weather station within the study area (Polo et al., 2010).The event is characterized by two main variables: D, duration (days) and P event , cumulative precipitation over the study area (mm).Different snow variables selected from the long term analysis results were also obtained for each event in the resulting 40-year event series at the region.The extreme data analysis was done on the subset events over the 95th percentile of P event .

Results
This section shows the different results of the long-term and extreme data analysis performed over selected annual and decadal variables, average over the study area (4584 km 2 ).The distribution over the five identified regions in the area (Fig. 1) is also presented to further asses some of the results.

Average annual precipitation and temperature during 1960-2000
On a regional basis, the 40-year mean annual values of P and T are 510 and 93 mm year −1 , for the total precipitation and snowfall, respectively, and 26, 12.5 and 0.4 • C for the average daily maximum, mean and minimum temperature, respectively.Figure 2a and b show the evolution during the reference period  of the annual precipitation (total precipitation and snowfall) and the maximum, mean and minimum daily temperature averaged over the study area by the model from the available daily datasets measured at the weather stations in Fig. 1.The decadal trend for these variables has also been included (Fig. 2c, d).
proc-iahs.net/369/157/2015/Proc.IAHS, 369, 157-162, 2015 Globally decreasing annual rates were obtained for both the total precipitation and snowfall (Fig. 2a), with values of 4.136 and 1.250 mm year −1 , respectively.Their decadal cumulative values analysis (Fig. 2c) shows the same trend, with snowfall dropping from 108 mm year −1 in the first decade to the final decadal 75 mm year −1 .A greater dispersion of both annual variables can be observed in this last decade.
In the case of temperature, an increasing trend of the annual mean of the daily mean and maximum values was found (Fig. 2c), with global rates of 0.035 and 0.021 • C year −1 , respectively.On the contrary, a decreasing global rate of 0.009 • C year −1 can be observed in the minimum values.The same trends were found in the decadal analysis (Fig. 2d), without any clear pattern in the dispersion of the annual values for each decade.

Snow-related variables
The evolution of SC mean and SC max at the study region is represented in Fig. 3a, which 40-year means values are 0.055 and 0.38 m 2 m −2 , corresponding with areas of 254.02 and 1736.50 km 2 , respectively.A decreasing trend is found for the mean value (0.0007 m 2 m −2 year −1 ) against the growing trend found in the case of the maximum value trend (0.0007 m 2 m −2 year −1 ). Figure 3b and c show the decadal analysis for both variables; where similar trends are also observed, with values over 0.35 and between 0.05 and 0.06 for the decadal SC max and SC min , respectively.Nevertheless, a high dispersion is found between decades.Figure 4a includes the evolution of the annual number of days with precipitation and with snow over the whole study area.Mean 40-year values of 173 and 40 days are found respectively.Both variables experiment a global decreasing trend, which represents a loss of 42 days with precipitation and 11 days with snow over the area during the analysed period.Figure 4b and c also show decreasing trend with variable dispersion between decades, which for example is very high in the last decades for days with precipitation.
The evolution of the annual mean of the maximum and mean daily SWE values is represented in Fig. 5a, in which a high variability can be observed, with mean 40-year values of 38.73 and 10.91 mm, respectively.The global decreasing trend obtained for both variables is higher for SWE max , which doubles the rate obtained for the SWE mean (0.625 against 0.235 mm year −1 ).The boxplots (Fig. 5b, c) with the decadal analysis also shows decreasing trends, with a high dispersion in both the first and last decades; this last decade shows the highest value for the SWE max in the reference period, and the widest interval of values for both SWE max and SWE mean variables.

Extreme data analysis
Following the extreme event definition in Sect.3, 82 extreme precipitation events were identified during the reference period, with a mean duration and cumulative precipitation of 12 days event −1 and 120 mm event −1 .Figure 6a shows the evolution of the annual number of extreme events throughout the study period, with a mean value of 2.1 events year −1 in a range of 0-5 events year −1 .The relation between the annual cumulative precipitation associated to extreme events and the annual cumulative precipitation for each year is included in Fig. 6b; again, a globally decreasing rate was found (0.006 mm mm −1 year −1 ).Finally, the torrential nature of these events can be analysed from Fig. 6c quantifying the rainfall intensity for the extreme event per year; a global increasing trend is found of 0.061 mm day −1 .

Discussion
The Mann-Kendall test performed to assess the significance of the trend shows that only the trends found for, P , snowfall, T mean , SC mean , n day with precipitation , SWE max and SWE mean are statistically significant (Table 1).
The results show a clear decrease of the annual snowfall on a regional basis (−1.25 mm year −1 ), which is also reflected in a reduction of the annual mean SWE (0.23 mm year −1 ).Therefore, the shift in the precipitation regime seems to be the most determining driver of the snow variables annual trends.The annual snowfall correlates linearly with P (r 2 = 0.50), but this correlation is very poor with T mean (r 2 = 0.18), although some increasing trend could be observed with decreasing annual mean daily maximum temperature values.This behaviour differs from the trends observed in the Alps, where precipitation and snow trends have opposing trends and seem to be more related with temperature regime (Laternser and Schneebeli, 2003).The annual fraction of snowfall over the total precipitation in the study  The annual mean of SWE max experiments a higher decrease than the daily SWE mean value.Nevertheless, this is not accompanied by a great reduction in the annual mean of the maximum snow covered area.Thus, an average regional reduction in the thickness of the snowpack must be associated to these peak moments, which supports the precipitation regime being the main responsible of the snow dynamics over this reference period.
The torrential events are more common in the last decades, as the intensity of the extreme events showed.This is also appreciable in the greater dispersion associated to the extreme proc-iahs.net/369/157/2015/Proc.IAHS, 369, 157-162, 2015 events variables found when compared to most of the variables during this last decade.The non-linearity of the snow behaviour makes it necessary a sound modelling of these regions if future scenarios assessment is to be performed.During this reference period, a global decrease in the annual precipitation rate and an increase in the annual mean temperature values has been quantified.The dataset of temperature and precipitation downscaling for the different IPCC scenarios over the study area for this reference period  overestimate both the precipitation and temperature annual rates of decrease and increase to −2.1 mm year −1 and 0.03 • C, and thus fail to reproduce the measured values.The direct use of these projected data sets is never advisable but for snow regions must be made with extreme caution.The combination of long term weather variables observations and snow modelling constitutes a sound alternative for historical time series trends over mountainous areas and is an adequate basis for correcting and simulating future scenario projections.

Conclusions
The joint analysis of weather and snow variables showed a decrease trend in the extent and persistence of the snow covered area over Sierra Nevada range.The precipitation rather than the temperature regime seems to be the most relevant driver on the snow regime forcing in Mediterranean areas over the reference period .This poses a constraint for rigorous scenario analysis in these regions, since the precipitation pattern is poorly approximated by climatic models in these regions and further assessment by means of a sound snow modelling is required.

Figure 1 .
Figure 1.Location of Sierra Nevada Mountain Range in Spain; limits of the protected areas: National Park (dark grey coloured area) and Natural Park (light grey coloured area); location of the different weather station (black dots); and limits of the five regions in which the study area has been divided: R1 -Adra, R2 -Andarax, R3 -Fardes, R4 -Genil and R5 -Guadalfeo.

Figure 2 .
Figure 2. Evolution of the annual cumulative precipitation (a) and the annual mean, maximum and minimum daily temperature (b), averaged over the study area.Boxplot of the decadal mean annual cumulative precipitation (c) and the decadal mean maximum, mean and minimum daily temperature (d); the central mark represents the median, the upper and lower edges of the box are the 25th and 75th percentiles, respectively, the whiskers extends to the most extreme data point not considered outlier, crosses are the outliers, and dots represent the decadal mean values.

Figure 3 .
Figure 3. Evolution of the annual mean and maximum daily fraction of snow covered area over the whole study area (a).Boxplot of the decadal maximum (b) and mean (c) daily fraction of snow covered area; the central mark represents the median, the upper and lower edges of the box are the 25th and 75th percentiles, respectively, the whiskers extends to the most extreme data point not considered outlier, crosses are the outliers, and dots represent the decadal mean values.

Figure 4 .
Figure 4. Evolution of the annual number of days with snow and with precipitation over the whole study area (a).Boxplot of the days with precipitation (b) and the days with snow (c); the central mark represents the median, the upper and lower edges of the box are the 25th and 75th percentiles, respectively, the whiskers extends to the most extreme data point not considered outlier, crosses are the outliers, and dots represent the decadal mean values.

Figure 5 .
Figure 5. Evolution of the annual mean of the maximum and mean daily SWE over the study area (a).Boxplot of the decadal maximum (b) and mean (c) daily SWE; the central mark represents the median, the upper and lower edges of the box are the 25th and 75th percentiles, respectively, the whiskers extends to the most extreme data point not considered outlier, crosses are the outliers, and dots represent the decadal mean values.

Figure 6 .
Figure 6.Annual evolution of the number of extreme event in the 40-year study period (a); relation between annual cumulative precipitation associated with extreme events and the annual cumulative total precipitation (b); and intensity for the extreme events (c).

Table 1 .
Results of Mann-Kendall test with a level of significance α = 0.05 realized for each one of the variables analyzed to evaluate the significance of the trends found.