Snow evolution in a semi-arid mountainous area combining snow modelling and Landsat spectral mixture analysis

This study proposes the use of both physically-distributed hydrological modelling in combination with satellite remote sensing images, to study the evolution of the snowpack in the Sierra Nevada mountains, in southern Spain. The snowmelt-accumulation module inside WiMMed (Watershed Integrated Management in Mediterranean Environment) hydrological model was employed, which includes the use of depletion curves to expand the energy and water balance equations over a grid representation. Snow maps obtained from spectral mixture analysis of Landsat images were used to evaluate this model at the study site. The results show a significant agreement between observed and simulated snow pixels in the area, which allows production of sequences of snow maps with greater resolution than the remote sensing images employed. However, some mismatches do appear at the boundaries of the snow area, mainly related to: (a) the great number of mixed pixels; and (b) the influence of the snow transport by wind.


INTRODUCTION
Water from the melting of the snowpack is a basic resource for human consumption, agriculture irrigation or hydroelectric production in mountainous areas.Snow dynamics is conditioned on one hand by snow accumulation, redistribution and ablation, and the associated meteorological driving processes, and on the other hand by the influence of rough topography, with strong gradients of exposure, vegetation and shadowing.All these factors can be highly variable in Mediterranean environments, which increases the heterogeneity of snow distribution and makes it difficult to monitor and measure its evolution.Furthermore, spatially distributed hydrological modelling has become one of the most commonplace techniques for its study.Besides reproducing the physical processes involved in snow dynamics, models can produce a spatial distribution of snow at the scales significant for the reproduced processes, i.e. small catchment at high spatial resolution, or large catchment at very small resolution.A great effort has been dedicated to develop physical snow models, with different examples in the literature (Jordan, 1991;Tarboton et al., 1994;Marks et al., 1999).In semi-arid regions, special considerations must be included (Herrero et al., 2009), such as highlighted from the analysis of a point snowmelt-accumulation model in Sierra Nevada (southern Spain).However, the extension of physical models to large areas in such regions requires the scaling-up from point to cell calculations; this is usually made through the use of depletion curves, which also captures subgrid variability (Luce et al., 1999), and the availability of distributed snow measurements, which is achieved by means of satellite-based snow observations, used in the calibration and validation steps (Schmugge et al., 2002).The snow variable most extensively employed in these processes is the snow cover fraction (SCF), due to the feasibility of obtaining it from remote sensing sources (Parajka and Blöschl, 2008).In semi-arid regions, the extremely changing conditions favour the particularly distributed pattern of the snow, which usually appears as medium to small sized patches.This is a limiting factor which is why Landsat imagery (30 × 30 m spatial resolution) is the most recommended for studying snow evolution over these areas (Marks and Winstral, 1999;Pimentel et al., 2012).
Obtaining snow maps from Landsat images is easy following the physical properties of the snow along the electromagnetic spectrum.The use of the Normalized Difference Snow Index (NDSI), which compares the visible and near infrared bands allows the discrimination between snow covered and non-covered pixels (Dozier, 1989;Herrero et al., 2011).However, when a subgrid classification is required, a spectral mixture analysis can be applied (Painter et al., 2009).These techniques assume that the reflectance measured by the sensor from a given cell is a combination of the individual reflectances of the different kind of surfaces present within the endmembers.
This work analyses the spatial distribution regime of snow in the Sierra Nevada mountains (Spain) at the subgrid scale by means of application of the point snowmelt-accumulation model developed by Herrero el al. (2009), and obtaining a 10-year series of snow cover fraction maps from a spectral mixture model applied to Landsat imagery.

STUDY SITE AND AVAILABLE DATA
This study site was in the Sierra Nevada Mountains, southern Spain.They are a linear mountain range parallel to the coastline of the Mediterranean Sea, with altitudes from 1500 to 3500 m a.s.l.The typical mountain alpine climate is modified by its proximity to the sea, which generates semiarid and tropical conditions in the surrounding area, mainly on the southern hillside closer to the sea.The study area selected was the Guadalfeo River basin upstream of Rules Dam, with an extent of 1057.3 km 2 , which has been previously studied and hydrologically modelled.Different weather station networks are available in this area; in this work, hourly and daily datasets of precipitation, temperature, radiation (shortwave and longwave), pressure, wind velocity, and relative humidity from the Agroclimatic Information Network of Andalusia (RIA) and the Guadalfeo Project Network (Herrero et al., 2011) were employed as input to the snow model.Specific interpolation algorithms were used in this abrupt topography, including, in the case of temperature and precipitation, the consideration of residual conditioned by the altitudinal gradients (Herrero et al., 2007) and, for radiation, topographic effects are also considered (Aguilar et al., 2010).

METHODS
The accumulation-snowmelt cycles during the ten years from 2004 to 2013 were simulated at the study site using the physical and distributed model developed and calibrated by Herrero et al. (2009Herrero et al. ( , 2011)).The selected 67 Landsat TM and ETM+ images from the same period were processed and snow maps were derived by using a spectral mixture model.The model performance was tested against these results.This section describes both the snowmeltaccumulation model and the snow maps obtained.

Snow model
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).
The model assumes a uniform horizontal snow cover surface 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 equations ( 1) and ( 2): 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.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 uR, uE, uW and uM are the advective heat rate terms associated with each of the mass fluxes involved in equation ( 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, which are empirical functions that relate the fraction of cell area covered by snow with some other snow state variables, and quantify the decrease in snow cover within the cell to take into account the reducing of the area implied in the energy-mass balance.

Landsat imagery
This section describes firstly the different levels of pre-processing needed to transform row image digital numbers into reflectance values, and secondly the method employed to obtain snow cover maps from these datasets.
Preprocessing is composed of four different, sequential steps: radiometric, atmospheric, saturation and topographic corrections.First, radiometric calibration by means of rescaling factors (Chander et al., 2009) to transform calibrated Digital Number in the original images into radiance values at the sensor was applied.Secondly, the atmospheric Dark Object Subtraction (DOS), was selected; this algorithm assumes that all atmospheric effects are represented by a blackbody in the image, which by definition must absorb all the solar radiation (Chavez, 1996), and together with this assumption, several simplifications on the reflectance physics equation, such as a Lambertian surface and a cloudless atmosphere, are done.After the atmospheric correction, the reflectance values obtained require additional corrections due to the saturation problems that usually appear in these snow areas associated with the radiometric configuration of the satellite.Based on the hypothesis of a high correlation between spectral bands for snow, a multivariate correlation analysis between bands was made to recover the snow saturated pixels values (Karnieli et al., 2004).Finally, a C-correction with the land cover separation algorithm (Hantson and Chuvieco, 2011) was employed to homogenize direct solar and non-solar illuminated areas in this rough mountainous terrain.
Snow maps Two methodologies were employed to obtain snow cover maps from Landsat imagery.On one hand, binary maps (cover/non-cover) were obtained based on the physical properties of the snow along the electromagnetic spectrum: snow has very marked different extreme values in the visible (high reflectance) and infrared regions (low reflectance).The use of Normalized Difference Snow Index (NDSI), which compares electromagnetic answers over these spectral regions in combination with two fixed thresholds: NDSI > 0.15 and TM1 > 0.06, allows the discrimination between both, covered and non-covered snow pixels.
where TMi correspond to the reflectance value in Landsat band i.On the other hand, linear spectral mixture analysis was employed to estimate the snow cover fraction within each pixel in all the images (Painter et al., 2009).This analysis assumes that the reflectance measured at the sensor, RS,λ , is a linear combination of the reflectance values from each spectral endmember, that is, pure surface covers present in the area: where Fi is the fraction of the endmember i; Rλ,i is the reflectance obtained after the distinct preprocessing steps of endmember i at wavelength λ; N is the number of spectral endmembers taken into account; and ελ is the residual error at λ associated with the fit of the N endmembers.The least-squares fit to Fi is solved by means of a reflective Newton method.
Three significant endmembers were identified in the study area, snow, vegetation (brush creeping vegetation) and rocks (phillite), their spectrum being obtained from a digital spectral library (http://speclab.cr.usgs.gov/).

Snow cover fraction maps
Snow maps during the 10-year series of Landsat images were obtained from both methodologies.Figure 2(a) shows the evolution of the global snow area at the study site throughout the period 2004-2013.Parallel temporal trends can be observed, with overestimated values of SCF when binary maps are employed, as expected, since there is a threshold over which this analysis results in a covered classification.This is particularly relevant during accumulation phases, with maximum differences usually obtained for the peak SCF values over a complete cycle.Melting stages generally exhibit closer values, especially during the final melting cycle of each hydrological year.A root mean squared error (RMSE) of 15.6 km 2 was obtained by comparing both kinds of results over the study period, which represents a notable improvement with using the spectral mixture analysis given that the maximum snow cover area found ranged from 370 to 130 km 2 .The snow distribution is shown in Fig. 2(b) for three selected dates during the study period, associated with different states throughout the snow season (beginning of accumulation, consolidated snow and beginning of melting).Moreover, Table 1 includes the fraction of area associated with cover (SCF = 1), non-covered (SCF = 0) and mixed (0 < SCF < 1) situations for these dates at the study site, for both classification algorithms.
Obviously, the covered/non-covered analysis does not discriminate mixed cells.The results show that the non-covered area is adequately represented by this algorithm; however, a significant portion of cells considered as snow-covered do not correspond to SCF = 1, with an associated area that can be close to a 40% of the estimated snow area by the covered/non-covered analysis.These mismatched pixels are usually located within border sites, as can be seen in Fig. 2(b), where the interaction between the melting snow and rocks is stronger.

Snow cover simulations
The 10-year study period was simulated by the distributed snow model previously described, with the calibration values obtained by Herrero et al. (2009).Figure 3 shows the comparison of the average SCF value in the watershed simulated over the study period and the average SCF value obtained from each Landsat image.A good level of accuracy of the simulation can be generally observed.A global RMSE of 0.05 m 2 m -2 , associated with 50 km 2 of snow area values, was found.
Figure 4 shows a comparison between the spatial simulated and observed distribution of snow throughout the watershed for the three selected images, including snow maps.As can be observed, again the border cells are not always correctly represented.Moreover, some snow free areas do appear inside the completely covered area that are not represented by the model; this characteristic snow distribution can be usually found and is generally due to wind effects.The model does not consider redistribution of snow by wind transport and therefore these mismatches are likely to appear when windy conditions between two consecutive images are not negligible.

CONCLUSIONS
The distribution of snow cover under highly variable conditions was estimated by spectral mixture analysis of Landsat images, which has proved to be a powerful tool for providing time map series to be used in snow modelling.The spectral mixture analysis allows an adequate estimation of the coverage of snow in each pixel of the area, which resulted in a global reduction of the snow cover area of close to 40% in the mixed identified cells when compared to a simple covered/non-covered classification analysis, this being especially important in these semi-arid mountainous regions where the particular snow dynamics favour the appearance of a great number of mixed pixels.
The calibrated model by Herrero et al. (2009Herrero et al. ( , 2011) ) has been validated by the good agreement between simulated and observed snow maps.The combination of these snow cover fraction maps with physical distributed snow modelling allows an indirect validation of other simulated snow variables, i.e. snow water equivalent and snow depth.The mismatching areas may correspond, among other effects, to locations and periods for which wind transport, not being considered in the snow model due to its complexity at high altitudes, are not negligible.

Fig. 1
Fig. 1 Location of Sierra Nevada Mountain Range in Spain, limits of the study area (black line), Rules Dam (white cross) and weather stations employed in the study (black crosses).

Fig. 2
Fig. 2 Comparison between the snow area resulting from the covered-non covered analysis and the spectral mixture analysis for the Landsat images used in the study area.a) Temporal evolution, b) Selected snow maps for 10 November 2011, 24 February 2013 and 10 May 2011 (up, covered-non covered classification, and down, spectral mixture analysis)

Fig. 3
Fig. 3 Evolution of mean SCF over the watershed.Black line represents simulated values and black dots are the estimates from spectral mixture analysis of Landsat imagery.

Table 1
Fraction values of each type of area considered in the analysis (completely covered, completely non-covered and mixed) obtained from both algorithms for the 10 Nov 2011, 24 Feb 2013 and 10 May 2011 Landsat TM images.