The Spatio-Temporal Dynamics of Water Resources (Rainfall and Snow) in the Sierra Nevada Mountain Range (Southern Spain)

Abstract

This paper describes the use of a unique spatio-temporally resolved precipitation and temperature dataset to assess the spatio-temporal dynamics of water resources over a period of almost seven decades across the Sierra Nevada mountain range, which is the most southern Alpine environment in Europe. The altitude and geographical location of this isolated alpine environment makes it a good detector of climate change. The data were generated by applying geostatistical co-kriging to significant instrumental precipitation and temperature (minimum, maximum and mean) datasets. The correlation between precipitation and altitude was not particularly high and the statistical analysis yielded some surprising results in the form of mean annual precipitation maps and yearly precipitation time series. These results confirm the importance of orographic precipitation in the Sierra Nevada mountain range and show a decrease in mean annual precipitation of 33 mm per decade. Seasonality, however, has remained constant throughout the period of the study. The results show that previous studies have overestimated the altitudinal precipitation gradient in the Sierra Nevada and reveal its complex spatial variability. In addition, the results show a clear correspondence between the mean annual precipitation and the NAO index and, to a much lesser extent, the WeMO index. With respect to temperature, there is a high correlation between minimum temperature and altitude (coefficient of correlation = −0.84) and between maximum temperature and altitude (coefficient of correlation = −0.9). Thus, our spatial temperature maps were very similar to topographic maps, but the temporal trend was complex, with negative (decreasing) and positive (increasing) trends. A dynamic model of snowfall can be obtained by using the degree-day methodology. These results should be considered when checking the local performance of climatological models.

1. Introduction

Mountains play an important role in the Earth System and the availability of water resources. They cover 25% of the global land surface, provide living space for 26% of the world’s population [1,2], supply freshwater resources to surrounding lowlands, and host a significant amount of biodiversity. Mountain environments are highly sensitive to changes in climate [3] and are sensors for early detection of climate change [4]. The high sensitivity of mountain environments to climate change has generated significant research to improve climate observations at high elevations [5]. The alpine regions of the Mediterranean are among the most climatically sensitive areas in the world [6]. In many areas of the world, climate monitoring in high altitude regions continues to be sparse and of low quality; observation networks are deficient and climate values are estimated by extrapolating from low-altitude data. This can be particularly misleading when extrapolating precipitation. There appears to be a general positive correlation between precipitation and altitude [7,8], but it is not as simple and direct as the high (negative) correlation between temperature and altitude [9,10].

The focus of the work presented here is a study on the spatio-temporal dynamics of precipitation and temperature in the high-altitude Sierra Nevada mountain range in Southern Spain in order to evaluate the spatio-temporal dynamics of water resources in the form of rainfall and snow. There have been many studies on precipitation in Southern Spain. For example, various aspects of precipitation in the Mediterranean region of Spain are considered in [11,12]. The authors of [13] studied daily precipitation in Spain for the period of 1951 to 2002 using data from 22 sites. The authors of [14] studied seasonal trends in precipitation in the Mediterranean Iberian peninsula from 1951 to 2000. The authors of [15] studied dry periods over the pluviometric gradient (from Gibraltar to Almeria) in Mediterranean southern Spain. The authors of [16] provided a database of gridded precipitation and temperature in Spain for the 1950–2003 period, but the grid spacing they used had a horizontal resolution of approximately 20 km.

These and other previous studies have estimated precipitation and temperature over large regions, such as the entire Mediterranean region or the whole of the Iberian Peninsula, without any altitudinal restriction. The work presented here focuses on precipitation and temperature within the high-altitude area of the Sierra Nevada mountain range, and specifically, on the connected region with altitudes from 1500 m above sea level (asl) to almost 3500 m asl, which comprises the most southern alpine environment in Europe. The objective was to build a database of daily precipitation and temperature within this high-altitude area and to interpret it by using various statistical analyses, as explained in the following sections. This database reveals the spatio-temporal evolution of water resources in the studied area.

2. Methodology

Estimating precipitation in mountainous areas is a problem of recent interest [17,18]. There are particular solutions to the general problem of evaluating the areal distribution of precipitation using limited precipitation records (see for example [19,20] for urban areas). However, geostatistics has proved to be the optimal technique for estimating daily precipitation and daily temperature using altitude as a secondary variable [21,22,23], among others. Altitude data are readily available from digital elevation models (DEMs) which provide complete coverage of an area. Rainfall data, however, are limited to sparsely located rain gauges. Temperature data are even more sparse, although they show a higher (negative) correlation with altitude. There are various geostatistical kriging procedures for incorporating altitude as a secondary variable. The most widely used is co-kriging [24], which is a multivariate geostatistical interpolator [25]. Essentially, the relative scarcity of the directly measured data for the primary variable (precipitation or temperature) is compensated for by using its spatial correlation with the more abundant secondary variable (altitude). The application of co-kriging requires estimates of the direct variogram (or covariance) of precipitation (or temperature), the direct variogram (or covariance) of altitude, and the cross-variogram (or cross-covariance) of precipitation and altitude (or temperature and altitude). These three estimated variograms are modelled by fitting permissible models, of which the most widely used is the linear co-regionalisation model [26]. Co-kriging is used to estimate daily precipitation and temperature on a regular grid, and seasonal and annual averages are estimated by summing the days contained in the chosen time period (or support in geostatistical terminology) of estimation (e.g., one month, one season, one year). As the number of rain gauges and temperature stations available for any particular day may be small, a climatological variogram is used [27,28]. The shape and range of the climatological variogram remains constant, while the nugget variance and partial sill are updated with daily precipitation and temperature statistics [28]. The mathematical basics of co-kriging are considered next.

Co-kriging is a geostatistical method for optimal multivariate spatial interpolation [26]. In geostatistics, spatial variable $Z\left(u\right)$ at spatial location $u$, representing, for example rainfall, minimum temperature, maximum temperature, etc., is modelled as a random variable. The set of all random variables $Z\left(u\right)$ in region $\chi$ of the space, $u\in \chi$, comprises a random function or random field $Z\left(u\right)$. With $\chi \subset {\Re }^d$ and $d=2$, the problem is two-dimensional, as is the case in the work presented here for precipitation or temperature and elevation. Thus, $Z\left(u\right)$ represents precipitation or temperature and is the variable of interest that is to be estimated by co-kriging.

It is assumed that $Z\left(u\right)$ is second-order stationary with constant spatial mean

$$\mathrm{E}\left\{Z\left(u\right)\right\}=m_Z,$$

(1)

and the two-point statistics, the covariance, and the variogram functions depend only on vector $h$:

$$C_Z\left(u,u+h\right)=C_Z\left(h\right)=\mathrm{E}\left\{Z\left(u\right)Z\left(u+h\right)\right\}-m_Z^2,$$

(2)

$${\gamma }_Z\left(h\right)=C_Z\left(0\right)-C_Z\left(h\right),$$

(3)

$$C_Z\left(0\right)={\sigma }_Z^2.$$

(4)

where $m_Z,{\sigma }_Z^2,{\gamma }_Z\left(h\right)$ and $C_Z\left(h\right)$ are, respectively, the mean, variance, variogram, and covariance of the random function $Z\left(u\right)$ and $\mathrm{E}\left\{.\right\}$ is the mathematical expectation operator.

In the simplest form of co-kriging, a variable of interest, or primary variable (e.g., temperature or precipitation), is estimated based on the experimental values of the variable and the experimental values of a secondary variable (e.g., altitude) that is correlated with the primary variable. The co-kriging estimator of the precipitation (or temperature) at any given geographical location $u_0=\left\{x_0,y_0\right\}$, where $\left\{x_0,y_0\right\}$ are the easting and northing coordinates respectively, can be expressed as:

$$Z^{ }\times \left(u_0\right)={\displaystyle \sum }_{i=1}^n{\lambda }_i^{0}Z\left(u_i\right)+{\displaystyle \sum }_{j=1}^m{\beta }_j^{0}Y\left(u_j\right).$$

(5)

where:

$Z\left(u\right)$ is the primary variable, precipitation or temperature, at location $u$, and $Y\left(u\right)$ is the secondary variable, altitude, at location $u$.

$n$ and $m$ are the number of values of variables $Z\left(u\right)$ and $Y\left(u\right)$, respectively, used in the estimation in Equation (5). Usually, these data are inside a neighbourhood centred on estimation location $u_0$.

The optimal weights for the linear estimation in Equation (5) are obtained by solving the corresponding co-kriging system; see for example, refs. [25,29]. If only the primary variable is used, ordinary co-kriging reduces to ordinary kriging of the primary variable. The same applies if there is no correlation between the primary and secondary variables.

For co-kriging, the direct variograms of the two variables and the cross-variogram (or direct covariances and cross-covariance) between the two variables must be estimated from the experimental data.

The unbiasedness of the co-kriging estimator in Equation (5) implies that the mean estimation error is zero:

$$\mathrm{E}\left\{Z^{}\times \left(u_0\right)-Z\left(u_0\right)\right\}=0.$$

(6)

This is achieved by including the following conditions in the co-kriging system:

$${\displaystyle \sum }_{i=1}^n{\lambda }_i^0=1,$$

(7)

and

$${\displaystyle \sum }_{k=1}^m{\beta }_k^0=0,$$

(8)

The variance of the estimation error can be written as:

$$\begin{gathered} \mathrm{Var}\left\{Z^{}\times \left(u_0\right)-Z\left(u_0\right)\right\}={\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^n{\lambda }_i^0{\lambda }_j^0{C}_Z\left(h_{ij}\right)+{\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{j=1}^m{\beta }_i^0{\beta }_j^0{C}_Y\left(h_{ij}\right) \\ +{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^m{\lambda }_i^0{\beta }_j^0{C}_{ZY}\left(h_{ij}\right)+{\displaystyle \sum }_{j=1}^m{\displaystyle \sum }_{i=1}^n{\beta }_j^0{\lambda }_i^0{C}_{YZ}\left(h_{ji}\right) \\ -2{\displaystyle \sum }_{i=1}^n{\lambda }_i^0{C}_Z\left(h_{i0}\right)-2{\displaystyle \sum }_{j=1}^m{\beta }_j^0{C}_Y\left(h_{j0}\right)+C_Z\left(h_{00}\right), \end{gathered}$$

(9)

The co-kriging system is obtained by minimising the estimated variance in Equation (1), subject to the unbiasedness conditions, which, in matrix form [30], is:

$$C\lambda =B,$$

(10)

with:

$$C=\left[\begin{array}{ccc}\begin{array}{ccc}{C}_Z\left(h_{11}\right)& \cdots & C_Z\left(h_{1n}\right)\\ ⋮& \ddots & ⋮\\ C_Z\left(h_{n1}\right)& \cdots & C_Z\left(h_{nn}\right)\end{array}& \begin{array}{ccc}{C}_{ZY}\left(h_{11}\right)& \cdots & C_{ZY}\left(h_{1m}\right)\\ ⋮& \ddots & ⋮\\ C_{ZY}\left(h_{n1}\right)& \cdots & C_{ZY}\left(h_{nm}\right)\end{array}& \begin{array}{cc}1& 0\\ ⋮& ⋮\\ 1& 0\end{array}\\ \begin{array}{ccc}{C}_{YZ}\left(h_{11}\right)& \cdots & C_{YZ}\left(h_{1n}\right)\\ ⋮& \ddots & ⋮\\ C_{YZ}\left(h_{m1}\right)& \cdots & C_{YZ}\left(h_{mn}\right)\end{array}& \begin{array}{ccc}{C}_Y\left(h_{11}\right)& \cdots & C_Y\left(h_{1m}\right)\\ ⋮& \ddots & ⋮\\ C_Y\left(h_{m1}\right)& \cdots & C_Y\left(h_{mm}\right)\end{array}& \begin{array}{cc}0& 1\\ ⋮& ⋮\\ 0& 1\end{array}\\ \begin{array}{ccc}1& \cdots & 1\\ 0& \cdots & 0\end{array}& \begin{array}{ccc}0& \cdots & 0\\ 1& \cdots & 1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\right]$$

(11)

$$\lambda =\left[\begin{array}{c}{\lambda }_1^0\\ ⋮\\ \begin{array}{c}{\lambda }_n^0\\ {\beta }_1^0\\ \begin{array}{c}⋮\\ {\beta }_m^0\\ \begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\end{array}\end{array}\end{array}\right],B=\left[\begin{array}{c}{C}_Z\left(h_{10}\right)\\ ⋮\\ \begin{array}{c}{C}_Z\left(h_{n0}\right)\\ C_{ZY}\left(h_{01}\right)\\ \begin{array}{c}⋮\\ C_{ZY}\left(h_{0m}\right)\\ \begin{array}{c}1\\ 0\end{array}\end{array}\end{array}\end{array}\right]$$

(12)

where ${\mu }_1$ and ${\mu }_2$ are Lagrange multipliers or parameters that are used to include the constraints given in Equations (7) and (8).

The solution of the co-kriging system:

$$\mathsf{\lambda }=C^{-1}B,$$

(13)

provides the weights required in the estimator in Equation (5).

In the geostatistical literature, the form of co-kriging summarised above is known as “ordinary co-kriging” to distinguish it from “simple co-kriging”, in which the mean $m_Z$ in Equation (1) is known [26,29,31].

3. Case Study

The isolated Sierra Nevada mountain range (Figure 1), in the Andalucía region of southern Spain, has an E–W orientation and dimensions of around 70 km in the E–W direction and a mean of 20 km in the N–S direction, although its width decreases from West to East. The total area of the mountain range is around 1360 km². The Mulhacen Peak, at 3479 m asl, is the highest altitude in the Iberian Peninsula and is located at longitude 3.312° W and latitude 37.053° N. This mountain range is at an ideal location to monitor climate change because it is an isolated alpine environment at, for Europe, a relatively low latitude.

Almost seven decades (67 years from 1951 to 2017) of precipitation and temperature data, measured at 247 stations (rain gauges and temperature stations), are shown in Figure 1B. Not all the stations were in operation for the entire period, and the network has changed over time. For example, Figure 2 shows the evolution over time of the mean number of operating precipitation (rain gauges) and temperature stations. The number of operating rain gauges increased steadily during the 1950s and 1960s and reached a mean number of around 70 rain gauges during the 1970s and 1980s, increasing thereafter to a mean of 90 rain gauges from the 1990s until the last year of the study, i.e., 2017. Over the last 20 years of this sequence, there was a high degree of variability in the number of operating rain gauges, but the total number was never less than 50. The number of operating temperature stations has also increased steadily since 1951, but usually their number is around a third of the number of rain gauges.

Daily precipitation and temperature were estimated for the 67 years of the study, which involved the spatial estimation of these variables for 24,471 days. The estimation network has a spacing of 500 m in each of the X and Y directions. The estimation network has 201 columns and 90 rows (Figure 1C) but the statistics were calculated for the area of the Sierra Nevada mountain range, restricted to the black polygon in Figure 1C and representing a connected area with altitudes higher than 1500 m. This altitude of 1500 m was selected on the basis that it delineates very well the isolated alpine massif of Sierra Nevada, as can be seen in Figure 1C.

4. Results

Geostatistical co-kriging was applied to the database of experimental precipitation and temperature data. The linear regionalisation model was used with a spherical variogram model with a range of 30 km and a nugget variance and partial sill that were updated each day according to the daily variance of precipitation and temperature, i.e., following the methodology of the climatological variogram. The result is a new database comprising daily precipitation and temperature (minimum, maximum and mean) for 67 years on a 500 m × 500 m spatial grid (Figure 1C).

Using the new database, various averages in space and time can be calculated to highlight the spatio-temporal dynamics of precipitation and temperature for the Sierra Nevada mountain range. As shown in Figure 1C, the mountain range has been defined by including all points of the estimation grid that have an altitude greater than 1500 m and that yield a connected polygon that clearly defines the Sierra Nevada mountain range.

The first important result is the mean annual precipitation in the Sierra Nevada for the period from 1951 to 2017, which is represented in Figure 3. Figure 3A shows a raster colour map and Figure 3B shows lines of equal precipitation (isohyets) overlaying the raster DEM in Figure 1C. The most surprising result is the strong orographic influence in the spatial distribution of the mean annual precipitation. The most striking fact is that the dome of maximum precipitation has an NW–SE orientation rather than the E–W main topographical orientation of the Sierra Nevada orographic dome, as reported in most previous studies. The second interesting result is that the maximum coincides with the summit of the range but is in the upper part of the main valley of the southern slope. The third surprising result is that there is a second relative maximum (2 in Figure 3A) in the eastern part of the Sierra, also located on the southern slope. There is also a second maximum of annual rainfall (3 in Figure 3A), but it is located outside, and northwest of, the Sierra Nevada mountain range.

The estimated map of mean annual precipitation (Figure 3A) is consistent with Figure 2 of Sumner et al. (2003) [32], which represents the average annual precipitation (mm) for 1964–1993, but not with the detail given in the work presented here. In addition, the Confederación Hidrográfica del Guadalquivir (Water authority for the management of the Guadalquivir river basin) reported a mean precipitation of 618 mm in the river basin of the Canales dam (with a surface of 176 km²). A similar estimate of 633 mm was obtained for the areal value for the same river basin in the map in Figure 3A.

Another interesting result is the map of mean precipitation for the 1951–2017 period, but instead of integrating time for the year, it is integrated for each season: Spring (April, May, June), Summer (July, August, September), Autumn (October, November, December), and Winter (January, February, March). The results are shown in Figure 4 and can be interpreted as a decomposition of the mean annual precipitation in Figure 3 to produce the mean precipitation by season. This map should be assessed in conjunction with the results shown in

Table 1. The mean annual precipitation in Sierra Nevada for the 1951–2017 period is 575 mm, with a standard deviation of 191 mm. The table shows the way in which this annual precipitation is distributed by months and by seasons.

Month	Mean Precipitation (mm)	% of Annual Precipitation	Season	Mean Precipitation (mm)	% of Annual Precipitation
January	68.0	11.8	Winter (JFM)	197.9	34.4
February	67.5	11.7
March	62.4	10.8
April	55.8	9.7	Spring (AMJ)	115.5	20.1
May	41.7	7.2
June	18.0	3.1
July	3.9	0.7	Summer (JAS)	39.1	6.8
August	7.0	1.2
September	28.1	4.9
October	58.7	10.2	Autumn (OND)	222.9	38.7
November	77.7	13.5
December	86.5	15.0

, which shows the mean precipitation by month and season. Table 1 shows that, climatologically, the driest month is July, with less than 1% of the annual precipitation, while the wettest month is December, with 15% of the yearly precipitation. August is also a dry month, with 1.2% of the annual precipitation, while June and September, with 3.1% and 4.9% respectively of the annual precipitation, can be considered intermediate months. The remaining months, from October to May, can be considered humid months that define the rainy season. The seasons, from wettest to driest, are Autumn, Winter, Spring and Summer, with 39%, 34%, 20% and 7% respectively of the annual precipitation. The mean annual precipitation for Sierra Nevada for the 1951–2017 period is 575 mm with a standard deviation of 191 mm; this implies a high variability from year to year.

Figure 4 shows that the wettest seasons (Autumn in Figure 4C and Winter in Figure 4D) reproduce the annual mean precipitation shown in Figure 3. The general spatial trend of the annual precipitation in Figure 3 is also reproduced in the other two seasons (Spring and Summer) but with the difference that the maximum precipitation in Spring occurred outside the target area in the north-western sector (number 3 in Figure 3A), while the maximum precipitation in Summer occurred in the second maximum of precipitation inside the target area, located in the eastern part of the Sierra (number 2 in Figure 3A). This behaviour was due to the orographic effect and the provenance of the humid winds, as explained in the discussion section.

Another important result is the mean annual temperature (minimum, maximum, and mean) in the Sierra Nevada for the period from 1951 to 2017, which is represented in Figure 5A–C, respectively. Because of the high negative correlation between temperature and altitude, the temperature maps in Figure 5 clearly resemble the map of altitudes in Figure 1C. The spatial distribution of temperatures can be calculated for any temporal interval as, for example, in Figure 5D, which shows the minimum January temperature for the 1951–2017 period.

Figure 6A shows the spatial variability of the total annual snow obtained by assuming that the daily precipitation when the minimum daily temperature was below zero degrees Celsius would fall in the form of snow. Figure 6B shows that the percentage of total precipitation that fell in the form of snow increased with altitude and, within the Sierra Nevada mountain range, varied from 4.2% to 78.7%, with a mean of 24.9%.

Another interesting result is the time series of the areal average of the annual precipitation for Sierra Nevada (black polygon in Figure 1C) for the 1951–2017 period, i.e., a period of 67 years or almost seven decades. The estimated time series is shown in Figure 7A. The mean annual precipitation of 575 mm is represented as a horizontal, large-dashed line. The high variability of annual precipitation over the years is reflected in a time series with a standard deviation of 191 mm. Despite this high variability, the annual precipitation decreased linearly at a rate of 33 mm of precipitation per decade. The negative slope is statistically significant (i.e., differs from zero), with a confidence level of 95%. The fitted linear trend is shown as a solid red line in Figure 7A. This decrease can be seasonally decomposed as a decrease of 15 mm in Winter precipitation, 12 mm in Autumn precipitation, and 6 mm in Spring precipitation, while there was no significant change in Summer precipitation. In addition, the years with the lowest precipitation coincided with droughts in southern Spain or in the entire Iberian Peninsula and have been marked with blue arrows in Figure 7A. It is clear from the Figure that the number of droughts increased with time, and it is estimated that a mean precipitation in Sierra Nevada of less than 440 mm is indicative of a drought in the region. A different aspect is the evolution of seasonality, defined here as the percentage of annual precipitation comprising Autumn precipitation. The result is shown in Figure 7B, in which it is clear that seasonality has remained unchanged over the 1951–2017 period. Seasonality has been defined as the percentage of annual precipitation occurring in Autumn, because Autumn is the main wet season in the area, recording, on average, 40% of the precipitation of the whole year.

Although Figure 7A shows the time series of the mean total rainfall over the Sierra Nevada mountain range, the same time series can be calculated on a pixel basis to show the spatial variability of the slope of the total rainfall as in Figure 8A. It can be seen that inside the Sierra Nevada limits, the slope is always negative (rainfall decreasing with time), but it is statistically significant (at the 95% confidence level) in the western part and on the northern side of the eastern part.

A spectral analysis of the time series of mean annual precipitation may reveal hidden periodicities or cycles in precipitation. The maximum entropy method [33,34] was applied to the time series in Figure 7A and the estimated power spectrum is shown in Figure 9. There are important cycles lasting around 12 years, related to sunspot cycles. A cycle of 83 years is also related to solar activity but is less reliable because it is longer than the length of the time series of the available data. The cycles of 6.5, 3.4, and 2.7 years are related to the North Atlantic Oscillation (NAO) index.

The previous suggestion that there may be a relationship between annual precipitation and the annual NAO index has been investigated further. Both variables (annual precipitation and the annual NAO index) are plotted in Figure 10A, in which there is an obvious correspondence of wet years with a negative NAO index, whereas dry years are related to a positive NAO index. This negative relationship can be more clearly seen by comparing the polynomial trends fitted to both time series (thick solid line in Figure 10A). This clear relationship between precipitation in the Sierra Nevada and the NAO index is considered further in the discussion section, as is the issue of the altitudinal precipitation gradient.

With respect to temperature, Figure 11 shows, for the Sierra Nevada mountain range, the time series of minimum, maximum, and mean temperatures for the 1951–2017 period. It is clear from the figure that the slope is negative, which implies that temperature is decreasing at a rate of around 0.1 °C every 10 years, which seems paradoxical in a global warming context. However, it can be shown that all three slopes are statistically significant with a confidence level of 95%. Again, it is instructive to assess the time series of the different temperatures for the different pixels, i.e., the spatial variability of the trend across space, as shown in Figure 12. Figure 12A shows minimum temperature, Figure 12C shows maximum temperature, and Figure 12B,D shows the statistical significance of the different slopes. For both minimum and maximum temperatures, there are positive and negative trends, together with areas that are not statistically significant (at the 95% confidence level), areas that have positive slopes that are statistically significant, and areas that have negative slopes that are statistically significant. For both minimum and maximum temperatures, the statistically significant positive slopes (increasing temperature with time) are in the eastern part of the mountain range, while the statistically significant negative slopes (decreasing temperature) are in the western part of the mountain range.

Finally, the equivalent product for snow can be obtained by using a simple snow precipitation and melting model. The adopted model is precipitation in the form of snow when the minimum temperature is below zero degrees Celsius. The melting model used is the temperature index of the snow melt model [35], formulated as:

$$\begin{gathered} \mathrm{If} F_1\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)\le F_2 \to \mathrm{Snow}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)=\mathrm{Precipitation}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right) \\ \mathrm{If}{F}_3\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)\le F_4 \to \mathrm{Snow}\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)=\mathrm{Snow}\left(\mathrm{i},\mathrm{j},\mathrm{t}-1\right)-F_5\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)F_6 \\ {F}_1,F_3,F_5\in \left\{T_{min},T_{max},T_{med}\right\} \\ {F}_2,F_4,F_6\in ℝ \end{gathered}$$

where $\left\{T_{min}, T_{max}, T_{med}\right\}$ are the daily minimum, maximum, and mean temperature, respectively.

$\left(\mathrm{i},\mathrm{j},\mathrm{t}\right)$: is the cell with spatial location (i, j) for day (t).

$F_2, F_4, F_6$ are real numbers. The first two are temperature thresholds and the third is a multiplicative melting factor.

The parameters were selected using a simple calibration that gave the following: $F_1=T_{min}$, $F_2=0$, $F_3=T_{min}$, $F_4=1$, $F_5=T_{min}$ and $F_6=5$.

Figure 13 shows the experimental snow thickness data given in [36,37]. The thin line represents the point measurements of the changing thickness of snow over a four-year period at a point located in the Refugio de Poqueira (Sierra Nevada). The thick line represents the pixel measurements of snow from the dataset used in the work presented here and the previous simple model. The changing thickness results are similar but not identical because of the different supports (point and pixel) on which they are measured.

Figure 14 shows the changing area of snow cover in the Sierra Nevada mountain range over the 1951–2017 period. Days on which the entire mountain range was covered by snow are recorded as 100% snow cover. The detail of the snow-covered area is given in Figure 15. Figure 16 shows the thickness of snow in the cell that contains the ski resort town of Pradollano.

5. Discussion

Three important issues in the spatio-temporal dynamics of precipitation in the Sierra Nevada have been identified in this study: (1) the importance of orographic precipitation [38] in understanding the spatial variability of mean annual precipitation; (2) the decreasing trend in the mean annual precipitation of 33 mm per decade; and (3) the NAO influence on mean annual precipitation. Other minor issues that have been revealed are (4) the persistence of seasonality over the almost seven decades studied and (5) the complex issue of the altitudinal precipitation gradient.

The precipitation patterns over the Sierra Nevada mountain range for the period of this study can be interpreted physically by considering the orographic effect and the main provenance of humid flows from the Atlantic and, to a lesser extent from the Mediterranean, that affect the Iberian Peninsula in general and the study area in particular. The authors of [12] distinguished atmospheric patterns with a clear distinction between Atlantic and western Mediterranean disturbances that produce characteristic precipitation patterns over southern Spain. Most of the annual precipitation is drawn from the Atlantic Ocean and is mainly influenced by North Atlantic climatic processes [39], such as the polar cyclone (Iceland Low), and hence the North Atlantic Oscillation (NAO) [40]. Annual precipitation is sourced to a lesser extent from the Mediterranean and is mainly influenced by the Western Mediterranean Oscillation (WeMO) [41]. Isotopic studies [39] revealed two main sources of humidity that influence precipitation in the Mediterranean part of the Iberian Peninsula. Convective precipitation events from the Mediterranean are isotopically enriched and prevail during Summer (Figure 10). In contrast, precipitation with lower isotopic values is transported along Atlantic storm tracks, which dominate during Winter [42]. The authors of [12] distinguished atmospheric patterns (synoptic types) with a clear distinction between Atlantic and western Mediterranean disturbances that produce characteristic precipitation patterns in Mediterranean Spain. In addition, from the synoptic types they obtained, The authors of [12] summarised the main general scenarios that produce significant precipitation in Mediterranean Spain. The first scenario is a large-scale disturbance located to the west of the Iberian Peninsula and producing humid Atlantic flows that induce precipitation in western Andalucia. A second scenario is the passage of cold fronts over the Iberian Peninsula associated with higher latitude, low pressure systems that encourage precipitation in the mountainous areas of Andalucia. In a third scenario, relatively small lows at 500 hPa are found in the southern part of Spain, and the associated low-level flux over the Mediterranean from the east-southeast is warm and humid. This configuration leads to precipitation over the eastern flank of Spain [12]. These three scenarios generated the three maxima in mean annual precipitation shown in Figure 3A.

The observed decreasing trend, at a rate of 33 mm per decade, in mean annual precipitation in Sierra Nevada for the 1951–2017 period is consistent with findings by other authors. A general decreasing trend in precipitation in southern Spain from 1960 onwards was observed by [43]. The authors of [32,44], among others, have also identified this decreasing trend of precipitation in southern Spain. Here, we have shown this decrease in the mean areal value of the high Sierra Nevada mountain range. The authors of [45] studied the evolution over almost 100 years of the Azores high (the “centres of actions” for the weather in the Iberian Peninsula) and concluded that blocking anticyclones have become more prevalent over western Europe, which could explain the decrease in rainfall over the Iberian Peninsula in general and in the Sierra Nevada in particular.

There are clear correlations between the NAO index and the annual precipitation in the Sierra Nevada and between the NAO index and Autumn and Winter precipitation. This is because 73% of the annual precipitation occurs in Autumn and Winter. To a lesser extent, there is a correlation between the precipitation in Sierra Nevada and the WeMO index. The positive phase of the NAO reflects below-normal heights and pressures across the high latitudes of the North Atlantic and above-normal heights and pressures in the central North Atlantic [46]. Hence, there is a northward shift of the axis of maximum moisture transport [47], and thus, there is a decrease in precipitation over southern Europe in general (including Spain) and in the Sierra Nevada in particular.

Despite the variability of the total annual precipitation and its decreasing trend, there was variability, although there was no trend, in the seasonality for the observed 1951–2017 period. Here, seasonality has been defined as the percentage of annual precipitation comprising Winter precipitation. The conclusion is that, irrespective of the trend in annual precipitation, Winter precipitation will be the main seasonal precipitation contribution to the total annual precipitation.

Finally, there is the issue of the complex spatial gradient of mean annual precipitation with altitude. The altitudinal gradient between locations A and B can be defined as:

$$G_{AB}=\frac{R_A-R_B}{H_A-H_B},\mathrm{with}{H}_A>H_B$$

where $R_A$ is the mean annual precipitation at location A and $H_A$ is the altitude of location A. Given that the order of the locations is $H_A>H_B$, a positive gradient means that precipitation increases with altitude. It is clear from Figure 3B that the altitudinal gradient of precipitation for Sierra Nevada will be complex. Because of the orographic effects shown in Figure 3, in the best case (considering the windward side of the Sierra Nevada mountains), there will be an altitudinal gradient in the mean annual precipitation of 15 mm/100 m. There will be a gradient of 10 mm/100 m in the mean annual precipitation for the best location on the lee side of the Sierra Nevada mountains. However, in general, the altitudinal gradients will be smaller if any two arbitrary locations are considered. These altitudinal gradients in precipitation are significantly different to some of the values reported in the literature for the Sierra Nevada mountain range, e.g., 170 mm/100 m reported by [48], among others.

With respect to the spatio-temporal dynamics of temperature in the Sierra Nevada, there are several points that have been identified in this study:

(1)

The spatial distribution of temperature is highly correlated with altitude, with correlation coefficients of 0.85 and 0.9 for minimum and maximum temperature, respectively.

(2)

The annual temperature trend has a negative slope but, as it is not statistically significant, the annual temperature in the Sierra Nevada cannot be considered to be decreasing, even in terms of the minimum, maximum, or mean temperatures.

(3)

The spatial variability of the trend in annual temperatures is complex. There are areas that show an increasing trend while others show a decreasing trend and others are undefined, as the trend is not statistically significant.

(4)

By using a simple degree-day model, it is possible to evaluate the evolution of the Sierra Nevada snowpack by assessing a large number of spatial and temporal statistics. For example, the transient seasonal snowline could be evaluated.

6. Conclusions

A significant database of daily precipitation and temperature on a regular square grid with a spatial resolution of 500 m has been obtained by geostatistical co-kriging, together with an important instrumental record of precipitation and temperatures. The study area for this work is the Sierra Nevada mountain range in southern Spain, restricted to the area with altitude above 1500 m asl, which corresponds to what is generally understood to be the Sierra Nevada mountain range. The main results, based on the precipitation and temperature dynamics of the Sierra Nevada mountains and obtained by a statistical analysis of the high-resolution precipitation database, are the importance of orographic precipitation and the influence of the NAO on the total amount of annual precipitation. The orographic precipitation explains the spatial distribution of the mean annual precipitation with maximum precipitation on the windward side of the mountain. This implies that the highest precipitation is not at the summit of the mountain, as had previously been assumed and which has led to over-estimates of the altitudinal gradient of precipitation for Sierra Nevada. Seasonality has not changed over the 1951–2017 period, for which there are almost seven decades of precipitation records for the most southern alpine environment in Europe. This implies that, irrespective of the total amount of precipitation, most of it will fall in Autumn and Winter. Smith (2019) [49] anticipated that mountain meteorology research would require special attention because of the altered distribution of precipitation, orographic effects, and physical mechanisms by which mountains modify precipitation patterns in different climate zones. Temperature (minimum, maximum, and mean) shows a more predictable spatial behaviour because it is well correlated with altitude. However, the spatial variability of the temperature trend is complex, which makes it difficult to draw any conclusions about global warming, at least in this area. To address this issue would require a more exhaustive monitoring of mountainous areas and their water resources.

The work reported here is a contribution to the climatological study of alpine environments and should encourage further research in the meteorology and climatology of high mountain areas which store strategic water resources.