A Snow Cover Climatology for the Pyrenees from Modis Snow Products

The seasonal snow in the Pyrenees is critical for hydropower production, crop irrigation and tourism in France, Spain and Andorra. Complementary to in situ observations , satellite remote sensing is useful to monitor the effect of climate on the snow dynamics. The MODIS daily snow products (Terra/MOD10A1 and Aqua/MYD10A1) are widely used to generate snow cover climatologies, yet it is preferable to assess their accuracies prior to their use. Here, we use both in situ snow observations and remote sensing data to evaluate the MODIS snow products in the Pyrenees. First, we compare the MODIS products to in situ snow depth (SD) and snow water equivalent (SWE) measurements. We estimate the values of the SWE and SD best detection thresholds to 40 mm water equivalent (w.e.) and 150 mm, respectively , for both MOD10A1 and MYD10A1. κ coefficients are within 0.74 and 0.92 depending on the product and the variable for these thresholds. However, we also find a seasonal trend in the optimal SWE and SD thresholds, reflecting the hysteresis in the relationship between the depth of the snow-pack (or SWE) and its extent within a MODIS pixel. Then, a set of Landsat images is used to validate MOD10A1 and MYD10A1 for 157 dates between 2002 and 2010. The resulting accuracies are 97 % (κ = 0.85) for MOD10A1 and 96 % (κ = 0.81) for MYD10A1, which indicates a good agreement between both data sets. The effect of vegetation on the results is analyzed by filtering the forested areas using a land cover map. As expected, the accuracies decrease over the forests but the agreement remains acceptable (MOD10A1: 96 %, κ = 0.77; MYD10A1: 95 %, κ = 0.67). We conclude that MODIS snow products have a sufficient accuracy for hy-droclimate studies at the scale of the Pyrenees range. Using a gap-filling algorithm we generate a consistent snow cover climatology, which allows us to compute the mean monthly snow cover duration per elevation band and aspect classes. There is snow on the ground at least 50 % of the time above 1600 m between December and April. We finally analyze the snow patterns for the atypical winter 2011–2012. Snow cover duration anomalies reveal a deficient snowpack on the Span-ish side of the Pyrenees, which seems to have caused a drop in the national hydropower production.


Introduction
The Pyrenees mountain range is located in southwest Europe at the northern edge of the Iberian peninsula (43 • N, maximum elevation 3404 m a.s.l.; Fig. 1).Because of the large amount of precipitation it receives, the Pyrenees range is the water tower for a region covering northern Spain, Andorra and south France.The headwaters of three major rivers in southwest Europe, namely the Ebro, the Garonne and the Adour rivers are located in the Pyrenees mountains.These rivers and their tributaries provide critical water resources for various economic activities, including hydropower generation and crop production in the irrigated lowlands.
As most of the winter precipitation falls as snow in the Pyrenees, the snowmelt is an important contributor to the river flow and shapes the hydrographs of the Pyrenean rivers (Lopez-Moreno and Garcia-Ruiz, 2004;Bejarano et al., 2010).Spring snowmelt is extensively used in downstream areas for crop irrigation during the growing season.Snowmelt is also stored in the many reservoirs located on the Pyrenean foothills.The main functions of these reservoirs are to supply runoff for irrigation in summer and to produce hydropower.For example, in the Ebro Basin there are 299 operating dams, according to the governmental inventory (Ministerio de Agricultura, Alimentación y Medio Ambiente, 2011).The unofficial inventory from the Spanish association of dams and reservoirs (Sociedad Española de Presas y Embalses, 2008) indicates that at least one-fourth of the dams in the Ebro Basin are used for irrigation and twothirds for hydropower generation (other uses include river regulation, aquafarming, water supply to urban and industrial areas, etc.).In addition, many hydroelectric plants in the Pyrenees also work without the use of a dam (run-of-theriver hydroelectricity).As a result, 33 % of the hydroelectricity power plants in Spain are located in the Ebro Basin, most of them being connected to the Pyrenean rivers (Fig. 2).The potential hydropower of the Ebro Basin represents 19.5 % of the total potential hydropower in Spain (Ministerio de Agricultura, Alimentación y Medio Ambiente, 2011).In France, the Pyrenean rivers are also highly exploited for irrigation and hydroelectricity (Fig. 2).The Garonne River is known as the only large watershed in France where a structural imbalance between water resources, the needs of different users and aquatic environments is officially recognized (Dupeyrat et al., 2008).The rising pressure on the water resources in the Ebro Basin is also an important concern (Milano et al., 2013).
Apart from these hydrological services, the snow cover is also critical for the tourism sector in the Pyrenees.In particular, ski resorts are an important source of income and local employment (Rived et al., 2013).Given the importance of the snow cover in the Pyrenees, it is necessary to monitor its evolution.The snow depth is recorded daily at 19 stations at least across the whole Pyrenees mountain range.Some stations were already installed in the 1980s but most of the data are available since the 2000s.Yet, these ground measurements are essential but insufficient to describe the snow cover dynamics in the various topographic and climatic contexts of the Pyrenees mountain range.On the other hand, remote sensing data are spatially consistent and therefore can be very useful to complement in situ measurements.
Since the early 1980s, it has been shown that space-borne sensors operating in the visible and near-infrared region of the electromagnetic spectrum are very effective to map snow cover (Bowley et al., 1981;Dozier and Marks, 1987;Baumgartner et al., 1987).As of today, only low-to mid-resolution sensors such as AVHRR, PROBA-V or MODIS allow for global observations of the snow cover at a daily time step (without cloud obscuration) with a spatial resolution of 1 km to 250 m.Higher-resolution snow cover maps (30 m) are typically extracted from the Landsat program images, but they provide data at a lower frequency (16 days), which is generally inappropriate to monitor the snow cover as large snow area variations can occur within a few days during the melt season (Rango, 1993;Gómez-Landesa and Rango, 2002).
A suite of snow products were derived from Aqua and Terra MODIS data and released in 2000 (Hall et al., 2002).The MODIS snow products are now widely used for hydroclimate applications in snow dominated regions.These products were generated using the SNOMAP algorithm, which primarily relies on the normalized difference snow index (NDSI).The normalized difference vegetation index (NDVI) is also used to improve snow detection in forested areas (Klein et al., 1998).There are differences between Aqua and Terra products as Aqua MODIS band 6 detectors are not functioning, while band 6 is used for the Terra product to compute the NDVI (Hall and Riggs, 2007).As a consequence, the NDVI test for forested areas is not activated for Aqua products.The successive versions of the MODIS snow products were compared with snow maps obtained from other sources at similar or lower resolution (i.e., relative validation, Hall et al., 2002;Klein and Barnett, 2003;Maurer et al., 2003;Simic et al., 2004;Rittger et al., 2013) and validated using ground snow measurements (i.e., absolute validation) (Klein and Barnett, 2003;Maurer et al., 2003;Simic et al., 2004;Ault et al., 2006;Parajka and Blöschl, 2006;Arsenault et al., 2014).Most of these studies were done using the Terra MODIS daily snow product (MOD10A1).More comprehensive reviews can be found in Hall and Riggs (2007) and Parajka et al. (2012).Despite the variety of methods used among these studies, the results led to the same conclusion that Terra MODIS snow cover products have a higher overall accuracy than snow maps derived from VEGETATION or AVHRR.The typical absolute accuracy of MOD10A1 is 93 % but depends on the land cover (Hall and Riggs, 2007;Arsenault et al., 2014).The accuracy is lower in forested areas (Simic et al., 2004;Parajka et al., 2012).Other important sources of misdetection are cloud/snow confusion (Rittger et al., 2013), the variation of the sensor viewing angle and the reprojection from the original swath data to the sinusoid grid (Arsenault et al., 2014).The accuracy of Aqua snow product is similar although less documented.A comparison with Terra snow maps indicated a lower accuracy at least in forested areas as expected (Hall and Riggs, 2007).Alternative snow cover detection algorithms were developed for MODIS reflectance data, which showed better accuracies (Sirguey et al., 2008;Rittger et al., 2013).However these codes are not operational globally and imply the processing of a large amount of data to get a snow map, which restricts the number of potential users.In contrast, MODIS snow products are available at the global scale through the Internet usually with a short delay of 8 days.
The comparison with ground data is a robust and precise method to validate satellite snow maps.However, a limitation is that the spatial representativeness of each measurement is generally unknown.It can be assumed that a ground observation is valid for the whole pixel above a certain snow depth threshold (Maurer et al., 2003).The other type of validation is performed using Landsat-derived snow maps as a reference, since Landsat high-resolution data allow for an accurate detection of snow even in mountains (Salomonson and Appel, 2004;Hall and Riggs, 2007).The advantage of this method is that there is no problem with the spatial representativeness of the reference data.Moreover, it allows evaluating MODIS snow maps in regions where the in situ station network for snow monitoring is insufficient.This technique was initially used to assess the SNOMAP algorithm (Hall et al., 1995(Hall et al., , 2002;;Klein et al., 1998).The comparison was limited to one or a few Landsat scenes with low or null cloud coverage.The small numbers of scenes allowed for a manual removal of clouds.Since 2009, the Landsat archive is freely accessible, which enables one to generate many snow maps over a large area for various periods of the year.This allows assessing MODIS snow products for varying snow conditions (e.g., fresh snow in winter, ripe snow in spring).Recently, Rittger et al. (2013) have taken advantage of this exceptional archive to validate the MODIS snow products.They found large errors during the snowmelt period and in forest areas and conclude that MODIS snow products still require to be used with caution.
The main objectives of this study are twofold.
-First, we aim to evaluate the value of the MOD10A1 snow product to generate a daily snow cover climatology in the Pyrenees for the period 2000-2013.
-Second, we use this new data set to characterize the variability of the snow cover in the Pyrenees.
The first objective is a necessary step even though many studies have already assessed the MODIS snow products accuracy in other mountainous regions.Indeed, continuing validation is important to make the information reliable, since the combination of topography, land cover and climate varies from one region to another (Rittger et al., 2013).In our case, the Pyrenees present a large physiographic variability in a rather small area and is under the influences of both Mediterranean and North Atlantic climates.Moreover, this assessment is important because the MODIS snow products are getting more and more attention from the regional agencies in charge of the water and tourism management in the Pyrenees (Parra et al., 2006).In particular, an important question for the water practitioners is the effective snow detection threshold, i.e., the value of the snow water equivalent or snow depth for which a pixel is statistically marked as snowcovered in a MODIS product.
We used both in situ and remote sensing data to assess both Aqua and Terra MODIS daily snow cover binary products (snow/no snow) in the Pyrenees.For the first time, we assembled a French-Spanish data set of continuous snow measurements from the Ebro Basin agency in Spain and Météo-France.This data set was used to validate the MODIS snow products and to determine the optimal snow detection threshold.Then, Landsat scenes over the Pyrenees corresponding to 157 dates between 2002 and 2010 were processed to generate an independent snow cover product.We did not focus on the discrepancies between the MODIS products and the Landsat or station data on specific dates or regions.We rather aimed at characterizing the range of uncertainties at the scale of the Pyrenees mountain range and across the snow season.
For the second objective, we implemented a gap-filling algorithm based on previous studies to generate a gap-free snow cover climatology from February 2000 to July 2013.This allows us to characterize the variability of the snow cover duration at the scale of the whole Pyrenees and its relationship with the topography.We show an application of this climatology to characterize the anomalous snow cover pat-terns during the 2011-2012 winter, which was particularly dry in the southern Pyrenees (San Ambrosio et al., 2013).

In situ data
We assembled two important data sets of in situ snowpack monitoring in the Pyrenees (Table 1).Météo-France provided the snow depth observations at eight stations in the French Pyrenees from January 2000 to December 2010.The snow depth was recorded daily at 06:00 UTC with a 1 cm resolution.The Ebro Basin agency (Confederación Hidrográfica del Ebro) provided the snow depth and snow water equivalent data from 11 telenivometers in the Spanish Pyrenees managed as part of the ERHIN program (Evaluación de los Recursos Hídricos procedentes de la Innivación; study of winter water resources) (Parra et al., 2006).Each telenivometer is equipped with an acoustic snow gauge and a cosmic ray detector for snow water equivalent sensing (Paquet and Laval, 2006).The snow season is generally shorter at Météo-France stations (e.g., Aulus-Les-Bains, Bagneresde-Luchon, St-Lary-Soul) because they are located at lower elevations than the telenivometers (Table 1).

Landsat
We used the data acquired by Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+).The data were collected from the United States Geological Survey (USGS) and the European Space Agency (ESA).The Landsat scenes spanning the Pyrenees are numbered 200-030 to 197-030 in the Worldwide Reference System 2 (eastward).There are 157 dates in our data set distributed between January 2002 and December 2010 for which at least one of these Landsat scenes is available.

MODIS snow products
MOD10A1 (Terra) and MYD10A1 (Aqua) snow products version 5 were downloaded from the National Snow and Ice Data Center (Hall et al., 2006(Hall et al., , 2007) ) for the period 1 September 2000-2 July 2013.This corresponds to 4688 dates among which 4625 dates are available for MOD10A1 (98.7 %) and 3996 dates for MYD10A1 (85.4 %) since Aqua was launched in May 2002 and operational in July 2002.From this archive, 157 MOD10A1 were available on the same day as the Landsat data set for MOD10A1 and 139 dates for MYD10A1.
We also downloaded the MOD10A2 product, which provides the maximum snow extent from MOD10A1 over a compositing period of 8 days on the same grid.

Land cover
We used the Corine land cover 2000 raster data version 15 that covers both France and Spain.It is considered as reference data for land cover mapping at the Europe scale (Bossard et al., 2000).The Corine land cover was used to produce a map of forested areas by aggregation of the broadleaved forest, coniferous forest, and mixed forest classes (63 % of the study area).

In situ data processing
Some snow depth or snow water equivalent values from the telenivometer data set were negative in summer, probably because of a drift in the sensors' calibration factors.These negative values were set to zero.Otherwise the data were not filtered or corrected.

Landsat processing
The processing of a large number of Landsat images is only feasible with an automatic cloud detection algorithm because the cloud mask is not provided with the Landsat data.Here we applied a cloud detection algorithm developed for highresolution multispectral images from Landsat, Venµs and Sentinel-2 (Hagolle et al., 2010).
The Landsat data were processed as follows.
-Orthorectification: Landsat data available from USGS are already orthorectified.Hence, ESA data were orthorectified using USGS data as a reference following the orthorectification methodology of the Centre National d'Etudes Spatiales (Baillarin et al., 2004).All images were projected to Lambert-93.The superposition errors were 0.2-0.8pixels for USGS data and 0.3-0.9 for ESA data depending on cloud coverage.
-Radiometric calibration for both data sets was made following Chander et al. (2009).The radiometric accuracy is not critical for this application.
-The scenes were assembled and resampled to a 240 m resolution.
-The cloud mask and cloud shadows mask were retrieved based on the detection of abrupt changes in the reflectance time series for every pixel (multitemporal algorithm; Hagolle et al., 2010).
-An external mask of water bodies was applied (SRTM water body data) -The snow cover was detected based on the NDSI and the reflectance in the green and SWIR (shortwave infrared) channels (Dozier, 1989).The NDSI was defined as where ρ green (resp.ρ SWIR ) is the top of atmosphere reflectance in the Landsat green channel (respectively, shortwave infrared at 1.6 µm).A pixel is flagged as snow if the three following conditions are fulfilled: a. NDSI > 0.4 b. ρ red > 0.12 c. ρ SWIR < 0.16.
These criteria were applied on reflectance corrected for a first-order slope effect (cosine correction; Meyer et al., 1993).If a group of adjacent pixels was detected as snow but entirely surrounded by cloud-covered pixels, then these pixels were flagged as cloud, otherwise thick cold clouds may be detected as snow.The mask generated with these criteria was dilated with a circular radius of three pixels to improve the detection on the snow region boundaries, which have a smaller snow thickness and generally do not fulfill all the three previous conditions.These steps and the red and SWIR reflectance thresholds were adjusted from the original formulation (Dozier, 1989) based on a visual inspection of the fullresolution images over the study area.

MODIS snow products processing
We extracted the "Snow Cover Daily Tile" field from MOD10A1 and MYD10A1, which includes the snow/no snow and cloud masks.Our zone intersects the MODIS sinusoid grid tiles h17v04 and h18v04.The grids were assembled and reprojected with the nearest-neighbor method in Lambert-93 over a region of interest covering the Pyrenees using the MODIS reprojection tool (extent given in Fig. 1).The different classes in the original product were merged in three classes: no-snow (no snow or lake), snow (snow or lake ice), no-data (clouds, missing data, no decision, saturated detector).The MODIS snow masks corresponding to the Landsat dates were also resampled to a 240 m resolution in order to match our Landsat snow mask spatial resolution.
In a second phase of this work we implemented a gapfilling algorithm to interpolate virtually all the missing values from 1 September 2000 to 1 July 2013.The algorithm was derived from Parajka and Blöschl (2008) and Gafurov and Bárdossy (2009).It works in four sequential steps.
1. Aqua/Terra combination: for every pixel if no-data was found in MOD10A1 then the value from MYD10A1 was taken.Otherwise, the value in MOD10A1 was kept.
The priority was given to MOD10A1 because we found that MYD10A1 is less accurate (see Sect. 3).
2. Adjacent spatial deduction: each no-data pixel was reclassified as snow (no-snow) if at least five of the eight adjacent pixels were classified as snow (no-snow).
3. Adjacent temporal deduction: a no-data pixel was reclassified as snow (no-snow) if the same pixel was clas- sified as snow (no-snow) in both the preceding and the antecedent grid.The preceding and antecedent grids were searched within a sliding temporal window, whose size was incremented if there was still a no-data value.This means that if a pixel was marked as no-data in grid n + 1 and/or in grid n − 1 (sliding window of 2 days), then the algorithm started again with a sliding window of 3 days (i.e., the test was done with grids n + 1 and n − 2, or grids n + 2 and grid n − 1).For this study we allowed the window size to be incremented up to 9 days in order to fill a long gap from 20 March 2002 to 28 March 2002.The longest data gap lasted 17 days in summer 2001 due to a power supply failure of the MODIS instrument on board Terra (15 June 2001 to 2 July 2001).This gap was not filled because this period is not critical for snow cover monitoring.
4. The sparse remaining no-data pixels were reclassified using a classification tree (Matlab Statistics Toolbox; Breiman et al., 1984).For each date, the snow and nosnow pixels were used to fit a classification tree on four predicting variables derived from the geographic position and the topography.The variables were pixel elevation, aspect, easting and northing (i.e., x and y coordinates in Lambert-93).The tree was used to predict the class of the no-data pixels for this date.
The gap-filled pixels were flagged in the final product to maintain a record of the original data.

Comparison
We computed the confusion matrices between the MODIS product snow/no-snow classification and both reference data sets (i.e., in situ data or Landsat data), which were considered as the truth.Based on these results we computed the overall accuracy, precision and the kappa coefficient (noted κ; Cohen, 1960).The overall accuracy (AC) is the proportion of the total number of predictions that were correct (i.e., snow or no-snow).The precision (P ) is the proportion of the predicted snow presences that were correct.In the analysis, we mainly used the κ coefficient because this statistic incorporates both information on agreement and disagreement between the MODIS products and the validation data (Klein and Barnett, 2003).For in situ data, the comparison was performed both for snow water equivalent (SWE) and snow depth (SD) measurements.
1.In situ data: the MODIS product snow presence/absence was extracted at each station.For each date, a pixel was considered to be correctly classified as snow if the snow depth value on the same day (or the SWE value) is higher than a threshold value noted SD 0 (or SWE 0 ).We tested 40 logarithmically spaced values between 0 and 8 m for SD 0 and 0 and 3 m w.e.(water equivalent) for SWE 0 .The same method was applied to the SWE data from the telenivometers, the SD data from the telenivometers alone and the SD data from all stations.
2. Landsat data: every pair of MODIS and Landsat snow masks obtained on the same day was compared on a pixel basis.The comparison was made only for the pixels which were not masked by the union of MODIS and Landsat cloud masks (i.e., where snow or no snow detection was possible for both data sets).For MOD10A1, a total of 14.7 × 10 6 pixels were compared   with the Landsat snow masks, among which 13.4 % were classified as snow in the Landsat data.
The gap-filled product was also evaluated by comparison with the in situ data.The comparison with Landsat was not performed because most of the no-data pixels are due to cloud cover, which also obstructed the Landsat image on the same day.However, we compared our gap-filled product with the 8-day composite product MOD10A2, which is often used in hydrometeorological studies.

In situ data vs. MODIS products
A first visual inspection of the time series suggests that there is a good agreement between in situ observations and the MODIS product, as shown here in the case of Ordiceto station only (Fig. 3).
From these data an optimal SWE detection threshold is found between 20 and 60 mm w.e.(top row in Fig. 4).For SD, the range of optimal values is narrower between 100 and 120 mm (middle and bottom rows in Fig. 4).The optimum is nearly identical for MYD10A1 but the agreement is a bit lower than for MOD10A1 as indicated by the κ curve.In what follows, we have set SWE 0 = 40 mm w.e. and SD 0 = 150 mm (Table 2).The resulting confusion matrices (Table 2) and statistical measures (Table 3) further indicate that an excellent agreement is found between the SWE data and MOD10A1 (κ = 0.95), but the agreement decreases if MYD10A1 is considered.Another result is that the classification accuracy is higher with the SWE variable than with the SD variable for the same stations (Table 2).The agreement between MODIS and in situ data significantly decreases when considering SD for all available stations (Table 3).In any case, however, the agreement remains acceptable since the lowest κ is 0.74.
The seasonal differences in the detection threshold were further analyzed based on the SD data from the telenivometers and MOD10A1 (Fig. 5).The results show that there is a trend in the optimal SD value (Pearson's correlation R = 0.90, p = 0.002).The maximum κ are found for smaller snow depth in the early snow season than in the late season.Similar results were obtained with the SWE, but the relationship was slightly less significant (R = 0.86, p = 0.007).

Landsat vs. MODIS products
The results show that 81.5 % of the Landsat snow-covered pixels are correctly classified in MOD10A1 (Table 4).For Here only the telenivometers data and the MOD10A1 product were considered.The snow depth value giving the maximum κ was extracted for each month and plotted against the month (inset).The linear fit was calculated with log-transformed values of the SD (Pearson's correlation R = 0.90, p = 0.002).
MYD10A1 the result is similar (81.6 %), but the number of pixels that were compared is lower and the fraction of snowcovered pixels in this sample was also lower (Sect.2.2).As a consequence, the performances are slightly better for MOD10A1 (Table 5).Figure 6 illustrates the higher dispersion of MYD10A1 with respect to Landsat snow cover area.It also shows that the snow coverage tended to be underestimated by MOD10A1.The confusion matrices also indicate that MODIS has a lower snow detection than Landsat since the false negative rate is close to 18 % (Table 4).Yet, the κ coefficients indicate a good agreement for both MODIS data sets (MOD10A1 κ = 0.85, MYD10A1 κ = 0.81, both statistically significant at the 1 % level).There is no evident dependency between the comparison results and the observation season (Fig. 6), which suggests that the MODIS snow detection is not significantly deteriorated by the snow properties (e.g., lower reflectance of ripe snow cover in late spring).For both MOD10A1 and MYD10A1, there is a lower agreement with Landsat when the comparison is made only over the forested areas (Table 5).However, the loss of accuracy is higher for MYD10A1.This is consistent with the previous studies (Sect.1).In particular, the proportion of correctly classified snow pixels drops by 9.5 % (from 81.6 to 72.1 %) in forested areas, whereas it drops by 4.1 % for MOD10A1 (from 81.5 to 77.4 %).

Gap filling
A virtually gap-free snow cover product with a daily time step from 1 September 2000 to 1 July 2013 was generated from MOD10A1 and MYD10A1.The gap-filling reduced the fraction of no-data pixels from 49 to 0.38 % (Fig. 7).The fraction of no-data decreased by a factor of about 10-20 % at every iteration from the previous stage, except for the spatial deduction step, which had a low effect, and the classification tree, which reduced drastically the no-data fraction by 90 %.The classification tree allowed for a complete removal of the no-data values for all the dates in which some no-data pixels remained after the previous steps.The 0.38 % missing values correspond to the 17-day gap in June 2001 for which the classification tree was not applicable.We found a substantial agreement between this gap-filled product and all in situ snow depth data.The κ coefficient (κ = 0.75, N = 27 684) is a bit lower than the one obtained with MOD10A1 (κ = 0.79, Table 3) but nearly equal to the one obtained with MYD10A1 (κ = 0.74).
The total snow cover area from the gap-filled product was also controlled using MOD10A2 (Fig. 8).As an example, we show the hydrological year 2005-2006 to illustrate how the gap-filling made substantial changes to the initial data and The snow coverage is defined as the fraction of the cloud-free area which is covered by snow in the study area (Fig. 1).For each date, the cloud-free area is defined as the intersection of the cloud-free areas in Landsat and MODIS products.how it revealed the shape of a typical snow depletion curve during the melt season.Similar results were obtained for the other years.

In situ data vs. MODIS products
For the validation with in situ data, it was useful to merge the Ebro Basin and Météo-France databases because it enabled expanding the range in station elevation and thus to obtain a more robust conclusion.κ coefficients are within 0.74 and 0.92 depending on the product and the variable, but the highest accuracy was obtained between MOD10A1 (Terra) and the SWE measurements.As expected, the agreement was a bit lower when considering MYD10A1 (Aqua).The inclusion of Météo-France measurements resulted in a significant decrease of the agreement (Table 3).These stations have a lower mean elevation (Table 1) than the Ebro Basin telenivometers; hence, the snow cover is more discontinuous or "patchy" in their vicinity.For example, three stations from the Météo-France data set are located in valley bottoms (near hydroelectricity plants); hence, their spatial representativeness at the scale of the MODIS product pixel (about 500 m) is limited.In spite of these variations we could identify consistent detection thresholds in SD and SWE (Table 4, Fig. 4).However, a clear seasonal trend was detected in the optimal SD and SWE detection threshold (see Fig. 5 for the case of the SD).This result is interesting because it reflects the hysteresis in the relationship between the amount of snow on the ground and its extent, which was often observed in alpine catchments (e.g., Magand et al., 2014).Small snow depths can cover large areas during the accumulation period.However, the spatial variability of the snow depth increases over the snow season due to ablation and redistribution processes.As a consequence, the minimum snow depth value to cover a MODIS pixel increases over the snow season.

Landsat vs. MODIS products
Regarding the comparison with Landsat, we also obtained a good agreement, but a manual screening of the comparison maps revealed that a large fraction of the misdetections was located along the snow cover edges, in agreement with previous studies in alpine and arctic regions (Déry et al., 2005;Rittger et al., 2013).An example is shown in Fig. 9.The snow commission in MODIS products can be due to the effect of forest obscuration (Parajka et al., 2012) because the lower boundary of the snow cover is often situated in forested areas in the Pyrenees.In this study we could detect a deleterious effect of the forests by comparing MODIS products with Landsat over forested areas.This method assumes that the Landsat snow detection was accurate enough to be considered as a ground truth.This assumption may not be always be valid, although the snow classification method for Landsat is well established.Indeed we could also observe snow commission errors in our Landsat snow maps data set along the snow cover edges.This may artificially increase the agreement with MODIS products.However, we consider that Landsat misclassifications are less frequent than MODIS.
Another possible cause for this error, which was not specifically investigated here, is the effect of the MODIS sensor view angle.This is illustrated by the image of 19 March 2009 in Fig. 9. On this day, the Pyrenees are on the edges of both input granules from the MOD10_L2 swath product (i.e., the input data used to generate the MOD10A1), where the MODIS instrument "bow-tie" effect is the most pronounced (Gómez-Landesa et al., 2004).This configuration causes a distortion of the gridded snow product MOD01A1 over the Pyrenees.The consequence is an increase of the false negative along the edges of the snow cover.
Lastly, another likely explanation is that the surface temperature screening in the snow mapping algorithm (Hall et al., 2001) is too strict so it eliminates true snow pixels in low elevation areas.The thermal threshold was discarded for the reprocessing of the collection 6 of MODIS snow products (Hall and Riggs, 2013); thus, we can expect some improve- ments with respect to this issue at least in the next MODIS snow product release.

Gap filling
The cloud obscuration is an important drawback of snow products generated from remote sensing instruments operating in the visible-infrared wavelengths.Here we had to interpolate an important fraction of the pixels to produce a consistent snow cover data set (about 50 %; Fig. 7).Similar cloud cover was reported other studies.After combining Aqua and Terra snow products only 5.3 % of the cloud pixels were converted to snow pixels.However, this represents about the half of the total snow pixels in the final product.The largest reduction in cloud obscuration is obtained through the temporal filter for up to 5 days (Fig. 7), in agreement with previous studies (Parajka and Blöschl, 2008;Hall et al., 2010;Gao et al., 2011).Beyond 5 days a higher uncertainty in the snow maps is expected but it is a necessary tradeoff to further reduce the cloud obscuration in the Pyrenees.We examined the cloud cover in the original product to evaluate the spatial uncertainty due to the gap filling.We used MOD10A1 to map the probability of cloud occurrence in the study area over 2000-2013 (Fig. 10, top panel) (Fig. 10, top panel).The cloud probability in the Pyrenees is more important in the northwest because the prevailing westerlies bring moist air from the North Atlantic into the continent, whereas the southeastern Pyrenees are more influenced by the Mediterranean climate, with a lower nebulosity.The cloud cover map also reflects the rain-shadow effect due to the orographic lifting of the air masses coming from the Atlantic by the Cantabrian Mountains and the Pyrenees in the west coast of the Iberian Peninsula.The seasonal variability of the cloud cover also reflects the influence of both Mediterranean and oceanic climates (Fig. 10, bottom panel).The cloud cover probability decreases in summer but remains substantial throughout the year.The highest cloud probability is found in April.Unfortunately this one of the months when the snow cover monitoring is the most useful because it corresponds to the beginning or the middle of the snowmelt season.This is an issue especially if the MODIS products are to be used under real-time conditions for river flow forecasting.
5 Application of the snow cover data set

Spatiotemporal influences on the mean snow cover duration
The new gap-free snow cover data set was used to compute the mean monthly snow cover duration (SCD), i.e., the number of snow days in the Pyrenees.We represented here the mean SCD per elevation band to characterize the climatological influence of the elevation on the snow cover dynamics (Fig. 11).It shows that the number of snow days increases strongly from the 800-1600 m a.s.l.band to the 1600-2400 m a.s.l.band between November and April.This is consistent with López-Moreno and Vicente Serrano ( 2007 season occurs between March and June, in agreement with e.g., Lopez-Moreno and Garcia-Ruiz ( 2004).There is snow on the ground at least 50 % of the time above 1600 m between December and April.We further analyzed the snow cover duration variability as a function of the slope's aspect for the area above 800 m a.s.l.(Fig. 12).It shows that the snow cover tends to persist longer on north-facing and east-facing slopes.This is consistent with the expected effect of the solar radiation on the snowpack energy balance.North-facing slopes receive less solar energy.West-and east-facing slopes are exposed to the same insolation but west-facing slopes receive solar radiation in the afternoon at the hottest time of day, which explains why the snow melts faster than on the east-facing slopes.If we further normalize the SCD with respect to the mean monthly SCD (not shown here), we see more clearly that the difference between east-and west-facing slopes increases over the snow season (from November to June).This is consistent with the previous explanation, because the effect of the solar radiation becomes more evident during the ablation season.This process is contributing to the hysteresis in the relationship between the SWE and the snow cover content, which was identified above (Sect.4.1).
For this analysis we used the same elevation and aspect grid as for the gap filling (Sect.2.2.3).Hence, a fraction of the snow cover climatology was constructed based on the same elevation and aspect data.However, this fraction is too small to significantly modify the results of this analysis (Fig. 7).
This analysis was restricted to the mean snow cover duration; however, interannual variability in the snow line eleva- tion can modify these general features (Krajčí et al., 2011), as illustrated in the next section.

Snow cover pattern in winter 2012
As above, we used the new gap-free snow cover data set to compute the mean snow cover duration between 1 January and 1 March for each pixel over the period 2000-2013 excluding 2012.The snow cover duration between 1 January 2012 and 1 March 2012 was computed separately (Fig. 13).Both maps reveal an anomalous snow cover distribution during the 2012 winter.Above-normal mean snow cover duration is observed in the northern part of the Pyrenees, contrasting with a strong snow deficit across the southern part.In situ SWE measurements in the Ebro catchment corroborate this analysis, as a large deficit was also recorded in most stations (not shown here, but visible in the case of Ordiceto station Fig. 3).
Given the importance of the Ebro Basin for the hydroelectricity in Spain (Fig. 2), we extracted monthly energy production data in Spain from the Spanish Statistical Office database (Instituto Nacional de Estadistica, Ministerio de Industria, Energía y Turismo, 2013).The hydropower pro-duction dropped in early 2012 in comparison with the mean value over 1995-2012 (Fig. 14).The annual production was also lower in spite of a recovery in April-May.At the national scale the total production was higher than the 1995-2012 average (all sources of energy included), which means that the energy demand was high.Hence, it is likely that the 2012 drop in hydropower production was caused by the 2012 winter drought in the Spanish Pyrenees, which we also observed in the snow cover data (Fig. 13).Further analysis is necessary to establish if and how the snow deficit contributed to this drop.The gross of the snowmelt generally occurs between April and July, but snow melting can be important throughout the winter in the lower elevation areas.It is also possible that hydroelectric dam managers reduced their hydropower production in anticipation to the coming deficit of snowmelt in winter.

Conclusions
We found an overall very good agreement between the MODIS Aqua/Terra products and two independent snow cover data sets generated from in situ (stations) and remote sensing observations (Landsat).Landsat data confirmed that the snow cover edges are prone to commission (snow not detected), particularly when the sensor view angle is large.Also, the uncertainties in the final snow product increases due to the interpolation of the cloud-covered pixels, particularly in the northwestern Pyrenees and during the winter and spring months.
In spite of these limitations the results of this study support the conclusion that the MODIS snow products provide valuable information on snow cover at the scale of the Pyrenees range.
Using all in situ data we could determine a statistically optimal detection threshold, i.e., the snow depth or snow water equivalent value from which it is very likely that a pixel is classified as snow-covered in MODIS products.We found that an acceptable SWE detection threshold is between 20 and 60 mm w.e. and a SD threshold between 100 and 120 mm for both MOD10A1 and MYD10A1.We recommend to consider these ranges of values to convert the snow depth simulated by a snowpack model into snow cover area at the MODIS resolution in the Pyrenees, e.g., for model validation or data assimilation.
The MODIS snow products should be used more carefully for hydrology because it is less accurate in the transition areas where the snowmelt is fast.However, they can provide meaningful insights for climatological studies provided that the missing values are interpolated.This is particularly revealing for a transboundary mountain range like the Pyrenees where hydroclimatic observations are collected by various agencies without a joint framework or a common data depository.We used this gap-filled snow cover product to compute the climatological snow cover duration in the Pyre- nees for the first time, to our knowledge.This information can now be used to study the spatiotemporal dynamics of the Pyrenean snow cover since 2000.Here, we were able to reveal the asymmetrical snow patterns during the 2012 winter (Fig. 13).A strong snow cover duration anomaly is evident in the Spanish Pyrenees, reflecting a precipitation deficit which may have caused a temporary drop in the hydropower production at the national scale.In order to further use MODIS data to improve hydropower prediction, the best solution may be the assimilation of MODIS data into a hydrological model.However, an issue is the cloud cover which can reduce significantly useful data in real-time conditions.

Figure 1 .
Figure 1.Shaded relief map of the study area and location of the snow monitoring stations used in this study (map projection Lambert-93).The rectangle is the area of interest used for remote sensing data (corner coordinates: upper left x = 320 km, y = 6250 km, lower right x = 680 km, y = 6100 km).

Figure 3 .
Figure 3.Time series of in situ data and MODIS data at Ordiceto station.The black dots are the daily snow depth (SD) and snow water equivalent (SWE) measurements.The color bars in the background indicate the snow presence (blue), absence (green) and no data (white) from MOD10A1.

Figure 4 .
Figure 4. Comparison of MODIS products with in situ data.Each graph shows the variation of the statistical agreement between both data sets with the snow detection threshold.Left: the detection threshold is in SWE and evaluated with measurements from the telenivometers.Right: the detection threshold is in SD evaluated with all available stations (telenivometer and Météo-France).AC: overall accuracy, TP: true positive, FP: false positive, TN: true negative, FP: false positive, P : precision, k: kappa coefficient.

Figure 5 .
Figure5.Evolution of the κ coefficient with the value of the snow depth detection threshold for each month of the snow season.Here only the telenivometers data and the MOD10A1 product were considered.The snow depth value giving the maximum κ was extracted for each month and plotted against the month (inset).The linear fit was calculated with log-transformed values of the SD (Pearson's correlation R = 0.90, p = 0.002).

Figure 6 .
Figure6.Snow coverage calculated from Landsat and MODIS snow products (Terra MOD10A1, Aqua MYD10A1) for 157 dates (139 dates for MYD10A1) distributed between 2002 and 2010 over the Pyrenees.The color indicates the month of the year.The snow coverage is defined as the fraction of the cloud-free area which is covered by snow in the study area (Fig.1).For each date, the cloud-free area is defined as the intersection of the cloud-free areas in Landsat and MODIS products.

Figure 7 .
Figure 7. Evolution of the number of pixels classified as no-data (e.g., clouds) during the gap-filling procedure.

Figure 8 .
Figure 8. Snow cover area evolution in the Pyrenees over the hydrological year 2005-2006 for three products: MOD10A1, the gapfilled product from MOD10A1 and MYD10A1 (this study) and MOD10A2 (8-day maximum snow extent).

Figure 9 .
Figure 9.Comparison between MOD10A1 and Landsat for two dates.Left column: MOD10A1 vs. Landsat classification (TP: true positive, FP: false positive, TN: true negative, FP: false positive).MOD10A1 false negative areas are mainly located along the snow cover edges.Middle column: Landsat and MOD10A1 classifications used for the comparison.Right column: location of the input swath granules (MOD10_L2 product) used to generate the MOD10A1 tiles (sinusoidal map projection).The MOD10A1 snow mask is distorted on 19 March 2009 because it was constructed from the border areas of two input granules where the bow-tie effect is most pronounced.
Figure 10.(a) Fraction of cloud-covered pixels in MOD10A1 over 2000-2013.(b) Mean monthly cloud coverage from MOD10A1 in the study area.

Figure 11 .
Figure 11.Mean snow cover duration in the Pyrenees over 2000-2013 for different elevation bands.The snow cover duration is the number of days with snow in our gap-filled product.It was computed with the same elevation bands as in Fig.1, except for 3200-3400 m a.s.l.(only four pixels at the MODIS resolution).Otherwise the fractional areas of each elevation band in the study domain are 53 % (10-800), 31 % (800-1600), 13 % (1600-2400), and 2.5 % (2400-3200).

Figure 12 .
Figure 12.Mean snow cover duration in the Pyrenees over 2000-2013 for the four main aspect classes (W: west-facing slopes, N: north-facing slopes, E: east-facing slopes, S: south-facing slopes).The snow cover duration was computed only for the area above 800 m a.s.l.

Figure 13 .
Figure 13.Illustration of the anomalous snow patterns in the Pyrenees during the 2012 winter.Top panel: mean snow cover duration (days) in January and February.Middle panel: snow cover duration in January and February 2012.Bottom panel: snow cover duration anomaly.

Figure 14 .
Figure14.Monthly energy production in Spain (units: ktep).The right y axis is for hydropower and the left y axis is for the total energy production.

Table 1 .
Description of the in situ snow data used in this study.Stations 1-11 are telenivometers from the ERHIN program for which daily data of snow depth and snow water equivalent were available.Stations 12-19 are Météo-France stations for which daily data of snow depth only were available.Longitude and latitude are given in decimal degrees (WGS84) and elevation in meters a.s.l.

Table 2 .
Confusion matrices between the MODIS snow products and in situ data.The numbers correspond to percentages with respect to the number of in situ measurements.The comparison was performed both for snow water equivalent (SWE) and snow depth (SD) measurements.The detection thresholds were set to SWE 0 = 40 mm w.e. for SWE and SD 0 = 150 mm for SD.These thresholds yielded the best κ coefficients.

Table 3 .
Statistical measures of the comparison between the MODIS product and in situ data, computed from the results presented in Table2(N: sample size, A: overall accuracy, κ: kappa coefficient).

Table 4 .
Confusion matrix obtained from the comparison of Landsat snow maps with Aqua/Terra MODIS snow maps.The numbers correspond to percentage with respect to Landsat data.The results are given for all the study area (all), and the area covered by forests only (forest).

Table 5 .
Statistical measures of the comparison between MODIS products and Landsat, computed from the results presented in Table 4 (N: sample size, A: overall accuracy, κ: kappa coefficient).The results for the forested areas are indicated in parentheses.