Cloud obstruction and snow cover in Alpine areas from MODIS products

Snow cover maps provide information of great practical interest for hydrologic purposes: when combined with point values of snow water equivalent (SWE), they enable estimation of the regional snow resource. In this context, Earth observation satellites are an interesting tool for evaluating large scale snow distribution and extension. MODIS (MODerate resolution Imaging Spectroradiometer on board Terra and Aqua satellites) daily Snow Covered Area product has been widely tested and proved to be appropriate for hydrologic applications. However, within a daily map the presence of cloud cover can hide the ground, thus obstructing snow detection. Here, we consider MODIS binary products for daily snow mapping over the Po River basin. Ten years (2003–2012) of MOD10A1 and MYD10A1 snow maps have been analysed and processed with the support of a 500 m resolution Digital Elevation Model (DEM). We first investigate the issue of cloud obstruction, highlighting its dependence on altitude and season. Snow maps seem to suffer the influence of overcast conditions mainly in mountain and during the melting period. Thus, cloud cover highly influences those areas where snow detection is regarded with more interest. In spring, the average percentages of area lying beneath clouds are in the order of 70 %, for altitudes over 1000 m a.s.l. Then, starting from previous studies, we propose a cloud removal procedure and we apply it to a wide area, characterized by high geomorphological heterogeneity such as the Po River basin. In conceiving the new procedure, our first target was to preserve the daily temporal resolution of the product. Regional snow and land lines were estimated for detecting snow cover dependence on elevation. In cases when there was not enough information on the same day within the cloud-free areas, we used temporal filters with the aim of reproducing the micro-cycles which characterize the transition altitudes, where snow does not stand continually over the entire winter. In the validation stage, the proposed procedure was compared against others, showing improvements in the performance for our case study. The accuracy is assessed by applying the procedure to clear-sky maps masked with additional cloud cover. The average value is higher than 95% considering 40 days chosen over all seasons. The procedure also has advantages in terms of input data and computational effort requirements.


Introduction
Climate change is expected to have a significant impact on water availability of snow-dominated regions (Beniston, 2003;Barnett et al., 2005).The discharges generated by mountain areas largely contribute to streamflows of many rivers that cross the most populated and economically developed regions in Europe.Several authors have predicted a shift of the hydrologic regime of the main European rivers, due to a changed behaviour of the mountain valleys.Changes are predicted for the Rhine (Shabalova et al., 2003;Lenderink et al., 2007;Hurkmans et al., 2010;Junghans et al., 2011), the Rhone (Etchevers et al., 2002;Boé et al., 2009) and the Danube (Hagemann et al., 2009).Focusing on the Po River system, Alpine areas cover 35 % of the basin and contribute on average for 53 % of the total discharge (Vanham, 2012).Understanding the effects of atmospheric forcings on snow cover duration and distribution at the basin scale is a starting point for modelling climate change impacts on the hydrologic regime.Also, these variables can be P. Da Ronco and C. De Michele: Cloud obstruction and snow cover in Alpine areas from MODIS products considered reliable indicators of climate trends (Haeberli and Beniston, 1998;Beniston, 2003), since snow cover depths and distribution are not only the results of diurnal values of temperature and precipitation, but mostly the product of the history of these two variables from the beginning of the accumulation season (Beniston, 1997).
In this perspective, remote sensing can be useful for reconstructing recent changes of snow cover extent, distribution and duration in wide regions such as northern Italy.Moreover, the use of snow cover maps for hydrological purposes is an effective tool, since the combination of largescale information with local estimations or measurements of snow features makes it possible to estimate the snow water equivalent (SWE) stored within a river basin (Molotch and Margulis, 2008;Bavera and De Michele, 2009).Moderate Resolution Imaging Spectroradiometers (MODIS) employed by Terra and Aqua satellites provide a Snow Covered Area product (SCA) with 500 m and daily resolutions, which consists of binary maps whereby snow is detected at the pixel scale.The accuracy of MODIS Snow Cover Products depends on region, season, snow condition and land cover type (Klein and Barnett, 2003;Maurer et al., 2003;Simic et al., 2004;Zhou et al., 2005;Tekeli et al., 2005;Ault et al., 2006;Parajka et al., 2006;Hall and Riggs, 2007;Liang et al., 2008).In Europe, Parajka et al. (2006) compared daily MODIS snow maps with in situ data of 754 climate stations over the whole of Austria, reporting an average classification accuracy of 95 % on cloud-free days.The wide and heterogeneous Austrian territory presents similarities with our case study.Accordingly, we expect a similar quality of MODIS SCA product for northern Italy.However, the primary limitation in using MODIS snow products is that no information on ground conditions is available in areas hidden by cloud.During a year, clouds may obscure most of the study area restricting the potential of using such snow cover images.For example Parajka et al. (2006) indicated that, on average, clouds obscured 63 % of Austria in daily snow maps from February 2000 to December 2005.The possibility of benefiting from a reliable product with daily temporal resolution is therefore conditioned by the ability to estimate the presence of snow in overcast conditions.With this aim, several procedures for MODIS products have been developed and tested in different regions (Parajka and Blöschl, 2008;Gafurov and Bárdossy, 2009;Wang et al., 2009;Hall et al., 2010;Parajka et al., 2010;Paudel and Andersen, 2011).Such methods are based on a spatio-temporal combination of MODIS data and they can generate cloudless images having accuracy comparable to that of the source product.Contrary to the case studies in Gafurov and Bárdossy (2009) (Kokcha basin: elevations range from 416 m a.s.l. up to 6383 m a.s.l., about 75 % of the basin area lies above 2000 m a.s.l.) and Paudel and Andersen (2011) (Trans-Himalayan region: 96 % of the area lies in the elevation zone above 3000 m a.s.l., with 43 % of the area above 5000 m a.s.l.), many European rivers drain basins which cover altitudes from the lowlands up to 4000 m a.s.l. in the Alpine region.The majority of the drainage areas correspond to elevations where snow dynamics are fast and several snow accumulation-melt cycles may happen over the year.This makes snow mapping under clouds a challenge.A spatio-temporal combination of MODIS images for cloud reduction was tested by Parajka and Blöschl (2008) over the whole of Austria, a territory similar to our case.Testing some backward filters (2, 5 and 7 days), they obtained a minimum accordance with ground observation of about 92 % in relation to a 7-day temporal filter.The latter was applied after merging Terra and Aqua images.However, temporal filters leave a percentage of cloud cover depending on the allowed temporal window from which the procedure can derive information to add to the cloudy day.Moreover, the accuracy of this method was dependent on the season (Parajka and Blöschl, 2008).During the winter snow processes in mountains are slower and in midsummer snow maps tend to a steady layout with snow over glaciers and snowfields.In spring and autumn ground conditions change more quickly, due to snowfalls and snowmelt which determine several land/snow micro-cycles on the transition altitudes.On the other hand, spatial filters exist which use information of the same daily images to map snow in areas hidden by clouds.Some of them estimate land cover pixel by pixel, filling the gaps with data found in the nearest neighbours (Parajka and Blöschl, 2008;Gafurov and Bárdossy, 2009;Tong et al., 2009).Others are based on the concept of snow transition elevation, looking for the altitude above which all pixels within a region, a basin or a slope can be supposed to be snow covered (Parajka et al., 2010).Paudel and Andersen (2011) pointed out that snow processes are highly affected by local conditions such as slope, aspect and land use, and the limit of spatial filters falls with the adoption of data recorded in other areas for estimating land cover of each individual cell.
Relying on recent proposals (Gafurov and Bárdossy, 2009;Parajka et al., 2010;Paudel and Andersen, 2011), we developed an innovative cloud removal procedure based on five steps.The aim was to integrate and improve methods proposed previously in order to maximize the performance of the procedure with reference to the Po River basin.The novelty stems from the investigation of the advantages provided by the combined use of different temporal filters, coupled with the regional snow line method (Parajka et al., 2010), in comparison to their individual application.At the same time we avoided the requirement of data from other satellites or information on land cover types: our solution involves only Terra and Aqua daily snow cover maps and a DEM of the domain.A second novel aspect concerns the validation methodology.The test was performed masking several clearsky maps with fictitious clouds and applying the procedure to such additional obstruction.In the first stage, 25 MOD10A1 and MYD10A1 SCA products were considered individually and covered by layouts of cloud cover extracted from other days, just as in Gafurov and Bárdossy (2009).In the second stage, more continuous cloudless maps were covered by real patterns of atmospheric disturbance, so as to simulate the hardest conditions met by the procedure over the year.
This work is structured as follows.First, we present the stepwise procedure and the validation methodology, potentially applicable to other remote sensing data.The description of the procedure is preceded by a concise introduction which explains the order of the steps.Then, we provide a brief description of the case study, focusing on the most interesting aspects for our targets.Here, we include some details concerning the source MODIS products and the data pre-processing.Finally, we show and discuss the results obtained both in terms of effectiveness and accuracy.Here, the performances are compared to those of a mere 7-day backward temporal filter to show that the increased complexity is justified by the higher reliability obtained.A detailed examination of the accuracy of each individual step is provided as well.

Methods
Methods are presented here since they are potentially applicable to many binary snow products.However, in this section we will explicitly refer to MOD10A1/SCA and MYD10A1/SCA, Snow Covered Area maps from the processing of Terra and Aqua MODIS data.

Investigation on cloud obstruction
The issue of cloud obstruction is studied to better understand how much it would restrict the fruition of remote sensing products with daily temporal resolution.We reclassified the pixels of each snow map in three classes (cloud C, snow S, land L), and let N c,d (Z > z T ), N s,d (Z > z T ), N l,d (Z > z T ) be the numbers of cloud, snow and land pixels included in the map of day d, for elevations greater than z T .For z T = z min , the basin is considered in its entirety.
represents the percentage of area hidden by clouds for Z > z T , observed at the day d.In addition to daily values of cloud cover, average monthly, seasonal and annual values of C d are used in this study.For example, the average annual value C a is obtained as where n d is the number of recorded days for that year.Average seasonal values of cloud cover can highlight systematic trends of cloud obstruction, shedding light on the most affected periods of the year.The cloud removal procedure will be engaged especially in those periods.We conceived cloud-obstruction maps to highlight the relationship between elevation and the frequency of cloud cover.They represent images displaying the number of cloud observations per pixel within the domain, computed on annual or seasonal basis.This tool gives practical information on the elevation ranges where the cloud removal procedure is mostly employed.

Cloud removal procedure
The procedure is composed of five steps, summarized in Fig. 1 using a flowchart.The sequence order has been designed so that moving from step 1 to step 5 the temporal window of the data involved for cloud removing is progressively increased from 1 day to the entire year.The basic assumption is that data closest in time are more reliable for snow mapping despite cloud obstruction, especially when they are recorded both before and after the cloudy day.This hypothesis loses its foundation when the nearest information is too distant.In this procedure, data nearest in time are used when they come from observations of the previous 6 days.Thus, the threshold is set to 1 week (6 days plus the day under examination).In unfortunate cases when clouds hide the same area for more than a week, other approaches must investigate snow presence.We address these situations by detecting a likely ground cover on the basis of a year cycle composed of a land and a snow season.
The first three steps, as well as their order, are similar to those included in Gafurov and Bárdossy (2009).
Step 1 combines Terra and Aqua images of the same day, taking advantage of the fact that cloud cover may vary significantly even within few hours.Step 2 is a conservative temporal filter (+2/−2 days), which adopts data available in the nearest days but only when such observations are found both before and after the day under examination.Step 3 uses the regional snow and land lines (Parajka et al., 2010) and is based on snow cover dependence on altitude.It follows step 2 since previous steps make available a higher amount of valuable data that can be used for detecting snow distribution over elevations.
Step 4 consists of a 6-day backward temporal filter.Thus, it involves a week of MODIS images.The introduction of step 4 is justified by the good performance provided by the 7-day backward temporal filter used in Parajka and Blöschl (2008) for Austria.However, when clouds cover the same pixel for more than a week, it seems more likely that ground conditions are those typical of the season than those observed 7 (or more) days before.Here, step 5 intervenes, detecting two seasons (snow and land period) as part of a yearly macro-cycle.This final analysis considers 1 calendar year of snow maps.

Step 1
The first step crosses Terra with Aqua information of the same day.The base map is the daily binary product (SCA) in MOD10A1 provided by MODIS on board the Terra satellite.The choice of MOD10A1/SCA as the main product is a result of several validations it has been subject to (Hall and Riggs, 2007).Moreover, MODIS sensor on Terra is dedicated to land surface, being designed and calibrated to this target.Aqua MYD10A1 images are involved in order to classify pixels obscured by clouds in the source Terra images.The choice of source map is important since Terra and Aqua products may present differences due to the different times of satellite overpasses.Aqua observes the considered area a few hours later, in the afternoon, when cloudiness has usually increased and ground condition may have changed.Having chosen MOD10A1/SCA as the starting product, MYD10A1/SCA intervenes only by adding information for pixels not already identified as snow or land in MOD10A1/SCA.Thus, Terra maps rule out any incongruity.
We name Map T (d, x, y, z) = S, L or C respectively a snow, land or cloud cell in the Terra snow product, located at (x, y, z) on day d.Similarly we name Map A (d, x, y, z) = S, L or C a snow, land or cloud pixel within Aqua map.The first step returns Map 1 (d) composed of snow, land, and cloud pixels ready to be processed further.Since we are looking for a daily map and this passage merges observations of the same day, it actually does not introduce any loss of accuracy.

Step 2
The second step estimates land cover by deriving information from the neighbouring days, but only when such conditions can be considered stationary over an appropriate temporal window (Gafurov and Bárdossy, 2009).For each cloud cell Map 1 (d, x, y, z) = C, the algorithm analyses its classification in the snow maps of the previous 2 days and of the next 2 days.Map 1 (d − 1), Map 1 (d − 2), Map 1 (d + 1) and Map 1 (d + 2) may contain new information.When agreement is found between the classifications within a maximum time interval of 3 days, such condition is used for classifying the pixel of Map 1 (d).The classification rules are summarized in and Map 1 (d + 1, x, y, z) = L the algorithm does not apply any correction to the cell, leaving it to the subsequent steps to determine whether there is snow.The basic assumption which justifies this filter lies in the concept that if both the previous and the next days have the same surface cover, such condition is kept in cloudy weather (Gafurov and Bárdossy, 2009).Even though snowmelt may occur on overcast days, direct shortwave radiation is shielded by clouds that reduce the solar energy available for the melting process.On the other hand, when the area is snow-free during the previous and the next remote observations, it might be snow covered in the meantime, but snow amount and duration would not be significant for hydrological purposes.

Step 3
The third step is based on the elevation dependence of snow cover.The aim is to find the right altitude below which we can assume that the cloud mask hides land pixels, and the elevation above which all pixels are likely snow.On the one hand, temperature is usually linked to altitude by a linear relationship depending on temperature lapse rate.Snow distribution in mountain areas is the result of snowpack conditions evolving during the whole winter season; exceptional situations such as temperature inversion may occur, but they are not expected to produce a significant impact on snow patterns.On the other hand, temperature is frequently used to determine both the state of precipitation and the melt rates (Schaefli et al., 2005(Schaefli et al., , 2014)).The widely used temperatureindex melt model (Hock, 2003;Ohmura, 2001) is based on the assumption that melt fluxes are linearly related to air temperature, which is seen as an integrated index of the total energy available for melting.
Step 3 involves a criterion similar to the regional snow and land lines approach proposed by Parajka et al. (2010).The regional snow line, µ s , is the mean altitude of snow-covered pixels within the basin.The regional land line, µ l , is defined similarly as the average altitude of land pixels, calculated over the cloud-free part of the domain.In this work, step 3 adopts even a subdivision of the area per aspect classes.Depending on its exposure, each cell is matched with one of the four classes: N (315 ).Thus, we allow classes to have their own snow and land line altitudes, elevations that may be affected by solar radiation amount.This measure can help in reproducing accurate snow profiles during the melting seasons, when snowfalls are rare and direct shortwave radiation rules the melting process (Hock, 1999).Shortwave radiation may be adopted as the main representative of global radiation, since the other components depend on it, directly or indirectly (Dubayah and Rich, 1995;Kumar et al., 1997).Once individuated µ s and µ l , step 3 turns into snow the cloud mask standing above the snow line and into land the part below the land line: We have avoided a discretization of the domain in subbasins or homogeneous landscape units (Paudel and Andersen, 2011), to make the procedure independent of additional data.Other proposals define such threshold altitudes adopting the highest land pixel and the lowest snow pixel (Gafurov and Bárdossy, 2009;Paudel and Andersen, 2011).Here, the limitation results from the export of a local assessment, linked to a single cell wherever located, as information on snow distribution over the whole study area.Within a densely urbanized basin, isolated snow-free areas may stand at very high altitude also in winter, thus shifting upwards the assessment of the snow line.Further, any error of MODIS Snomap algorithms, which would classify isolate land cells where snow is expected (e.g. a snow area lying in shadow not individuated by the snow mapping procedure) or vice versa (e.g. the misclassification of cirrus clouds as snow in the period between May and October), may prevent the effectiveness of step 3.
We have introduced some limitations to the intervention of the regional line method.First, we allowed its application only in cases when the basin is at least 50 % cloud free, in order to ensure a representative estimation of the actual snow and land average altitudes.Furthermore, the snow line is calculated only if observed snow pixels are at least 5 % of the observed land pixels -5 % is a representative percentage of the snow pixels found in maps during the late spring on some cloud-free days.Throughout the winter, the ratio between the number of pixels classified as snow and the number of snow-free pixels should be higher even in the presence of clouds.If this condition is not met, few snow areas are observed and their mean altitude cannot be considered representative of the actual snow pattern for the whole study area, even when the sky is sufficiently clear.For the same reason, we let the algorithm skip the assessment of the snow line in summer (from June to September), when the percentage of snow areas is very low.The latter measure is designed to avoid errors arising from the classification of some clouds as snow; this misclassification may become significant in assessing the snow line, because of the lower number of actually snow-covered cells.

Step 4
Until this point, the request of maintaining snow maps with daily temporal resolution is satisfied.Steps 1 and 3 bring down cloud obstruction using information from the same daily map, while step 2 uses data from close days but only when ground conditions can be assumed constant over time.Further passages must solve the transition band which lies between land and snow regional lines, and all the other cases when step 3 is not involved.
Step 4 uses information derived from a selected time window.Ground conditions are thus approximated as stationary over a chosen period.The more we expand this window, the higher is the probability that ground cover is changed within (Parajka and Blöschl, 2008).Here, we chose a 6-day backward temporal filter.A pixel obscured by clouds on day d is classified to snow or land depending on the last time it was observed in the previous 6 days (i.e. the most recent snow/land classification is used given that it falls in the previous 6 days).The threshold of 6 days ensures a temporal resolution of the output map higher than the MODIS 8-day snow product, which contains the maximum snow cover extent over an 8 day window.Looking at the results provided by a 5-day and a 7-day filter for Austria (Parajka and Blöschl, 2008), 6 days seems to be a good compromise to ensure a good level of accuracy while reducing cloud pixels to a very low percentage.Since errors generated by steps 2 and 3 will be used in step 4, the 6-day filter removes remaining clouds exploiting output maps by step 1.This ensures that the most recent classifications come from observations and not from estimates.Otherwise a sensitivity analysis would be required.

Step 5
The remaining clouds are processed by a fifth step, named "seasonal filter", which outlines a yearly macro-cyclecomposed of a land and a snow (or accumulation) season.It was developed similarly to the last step of Gafurov and Bárdossy (2009).It works at the pixel scale individuating a land season and a snow season per pixel.Each cloud cell is then classified as snow or land in accordance with the calendar day.Several features of the location (e.g. the exposure to solar radiation and the land use) intervene in the snowmelt process; for these reasons, a pattern designed specifically for each cell can be an advantage in predicting the individual behaviours.
The first target is the right definition and identification of an accumulation season and a land season per pixel.As accumulation season we define a time period when a pixel is expected to be snow covered (generally enclosed within winter and spring).Even though land observations are possible meanwhile, it is assumed that in the accumulation months the pixel is likely covered by snow when still masked by clouds after step 4. On the other hand, once the "land season" has been individuated, each cloud observation not processed by the previous steps will be changed to land.Step 5 uses maps released by step 1, since we prefer that these two time frames come from snow cover actually observed from one satellite or the other.The snow-start and land-start flags are individuated based on an improved approach described as follows: -Given a pixel, its accumulation-start flag is individuated as the first day when snow is observed, providing that it is followed over time by other n s snow observations not interrupted by land.This implies that the accumulation season starts at the first snowfall, but only if the pixel is still snow covered or hidden by clouds in the next days, until the (n s + 1)th snow observation is met.The aim is to avoid false identifications of the accumulation-start flag due to occasional snowfalls that may occur, even at low altitudes.The threshold number of snow days n s is set depending on elevation, using three elevation bands, with a mean value of n s = 2 for the transition elevations.Such threshold decreases with altitude, since at lower elevations it is likely that the first snow cover is due to an occasional event, while in the higher areas one or two snow observations uninterrupted by land seem an appropriate signal of the beginning of winter (a complete view of the thresholds is provided in Table 2).
-The same rule is adopted for individuating the landseason start flags, in spring or summer (Fig. 2).Here, the threshold value n l is set increasing with elevation.This measure prevents situations in which exceptional absence of snow in the uplands (e.g.due to anomalous high temperatures) is taken as indicator of the absolute end of the snow period.
The choice of involving more days after the first land/snow observation protects from possible misclassification of MODIS algorithms, which may confuse some types of cloud with snow, and also where and when snow is not expected.In particular, the classification of some cirrus clouds as snow is a well-known issue (Parajka et al., 2006).We specify that a snow accumulation period is designed only for those areas lying at altitude higher than z T = 600 m a.s.l.-the threshold altitude was chosen based on some studies on snow duration in the Alps (Beniston, 1997;Valt and Cianfarra, 2010).As a consequence, pixels not already classified and located lower than z T are converted to land regardless of the time of the year.At low elevations long-lasting snow covers are infrequent, while several short cycles of accumulation-melt may occur.In these situations snow depths remain relatively thin, thus our choice may lead to a slight underestimation of the snow amount which would not justify the delineation of each brief snow cycle.
At high elevations, snow cover develops according to an annual cycle composed of two macro-periods, evolving from the accumulation to the melt season.However, the life of snow cover is highly affected by altitude, which provides an indication of the persistence of low temperatures.The Italian Alps are mostly distributed on transition altitudes, where several land/snow micro-cycles occur, characterized by snow durations of a few days to several months.Even with the new expedient, the separation of the year into two seasons per pixel may lead to rough errors of snow overestimation at the beginning of the winter season and underestimation after first (but maybe not last) melting completes.That is why step 5 is preceded by a temporal filter in step 4.
The added value of considering more snow/land observations will be presented in the results.
Step 5 is validated individually, considering this novel setup in comparison with the simple "first observation" approach (Gafurov and Bárdossy, 2009).

Validation methodology
We check the accuracy of the procedure improving the methodology proposed by Gafurov and Bárdossy (2009), and also used in Paudel and Andersen (2011).Some cloud-free days, henceforth considered as "truth", are selected and covered by wide cloud masks.Such cloud cover is borrowed from dates when both Terra and Aqua products register overcast conditions.The use of observed layouts of cloud cover ensures that the artificial images assigned to the procedure correspond to configurations that may occur.In order to avoid the possible no-likelihood of artificially masked maps to time series of clouds that can occur (e.g.several days of continuous cloud cover) testing days were chosen within time frames of cloudy days alternated with clear-sky days.Then, the procedure is applied to the maps, in order to reclassify the additional clouds.Finally, the removed-cloud product is compared with the native clear-sky Terra map.
In this stage we are not interested in the reliability of the source MODIS products, which is a widely studied issue (Hall and Riggs, 2007).Hence, this validation strategy presents some advantages with respect to the commonly used evaluation based on in situ data.For instance, one of the more discussed issues is whether a point value can be considered representative of snow coverage within a 500 m cell (Wang et al., 2009).Consequently, appropriate snow depth thresholds for in situ data must be chosen.Snow depths of 4 cm (Wang et al., 2009), 1 cm (Maurer et al., 2003) and 2.54 cm (Simic et al., 2004;Tekeli et al., 2005) have been used in the literature.This influences the results of the validation (e.g.Parajka and Blöschl, 2008;Gao et al., 2010).Furthermore, a validation carried out through snow-depth values is highly influenced by the number of stations available (and their altitudinal distribution), as well as by the temporal continuity and reliability of data.In Italy, few stations exist at altitudes over 2500 m a.s.l., where snow stands until the late spring and summer.In contrast, by masking clear-sky images with wide clouds, the number of pixels over which the test is performed is certainly higher.

1 day additional cloud mask
Once the cloudy day has been chosen, clouds from Terra are overlapped to the clear-sky Terra maps and the same is done for Aqua.
Given a clear-sky image on day d, fictitious cloud cover introduces an additional number of cloud pixels, N c , within the domain, whose percentage is The N c is calculated as the number of cloud pixels in the masked image that were cloud free in the source map.The A d,T and A d,A refer respectively to the additional cloud cover on Terra and Aqua images.The N tot represents the total number of cells included in the basin.Thus the percentage of cloud cover in the masked Terra map is where N c,T is the number of cloud pixels in the clear-sky Terra product at day d.Likewise, C d,A indicates the percentage of cloud cover for Aqua images after the introduction of fictitious clouds.The gap-filling procedure reclassifies N c,T either to snow or land.Comparison with the original Terra products identifies three possible situations: overestimation errors (land to snow L → S), underestimation errors (S → L) and agreement (snow to snow S → S or land to land L → L).Let a and b be the numbers of overestimations and underestimations, and c and d the numbers of pixel correctly reclassified as snow or land.The overall degree of agreement D A is (Paudel and Andersen, 2011) The degree of disagreement is The latter is due to overestimations and underestimations The overall accuracy is achieved by the average of the different performances, i, weighted on the percentage of pixels under examination per testing day: The variance, σ 2 , is calculated as

Multi-day additional cloud mask
Here, we test the procedure to deal with the worst-case scenario proposing a novel validation strategy.The idea is to simulate cloud obstructions which affect several days.Some groups of continuous days characterized by very low percentages of cloud obstruction are chosen within the period 2003-2012.In parallel, atmospheric disturbances covering the same number of days are extracted from other Terra and Aqua imagery and overlapped.This generates a spatiotemporal additional mask borrowed from the same season of another year.Consequently, the gap-filling procedure loses several days of data, exactly when more information on snow cover are available.At the same time the likelihood of the artificial images to match natural dynamic conditions is ensured.The degree of accuracy is calculated per validation day as done for the 1-day masks.

Case study and data set
The Po is the largest Italian river and one of the most important fluvial systems in Europe from an economic point of view.It flows more than 650 km eastward across northern Italy, from a spring seeping at Pian del Re, Piemonte, through a wide delta into the Adriatic sea between Venice and Ravenna.Its drainage area covers approximately 74 × 10 3 square kilometres, of which about 71 × 10 3 are located in Italian territory, a quarter of the total national territory.The remaining area is located for the most part in Switzerland (Toce River basin) and a smaller part in France.Given the level of utilization of water resources for agricultural and industrial purposes, the Po Basin is a focal point of the Italian national economy.Through its tributaries, the river drains mountain regions that reach altitudes above 4700 m a.s.l. in the Aosta valley (Fig. 3a).The hypsographic curve for the Po River basin is shown in Fig. 3b, where the x-axis represents the surface areas (or relative surface area) which lie above and below a marked elevation while the y-axis represents elevation above sea level.The curve has been obtained from a DEM with 500 m spatial resolution.Figure 3b shows that more than 30 % of the area lies above 1000 m a.s.l., where snow dynamics are expected to have an important effect on the hydrologic cycle.The spatial and temporal evolution of the snow cover over the basin is thus one of the core issues for characterizing the hydrologic regime of the river, as well as for assessing its susceptibility to changing climate conditions.
After the pre-processing of the data set (Sect.MODIS capability in monitoring the annual fluctuations of snow cover extent, Fig. 4a shows a comparison between snow observed by the Aqua MODIS on 28 June 2005 and a digital map of Lombardy's glaciers.The latter is the result of processing and comparison of areal data taken from orthophotos (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007), and from glaciological registers (more information available at https://www.dati.lombardia.it/Territorio/Ghiacciai-e-Nevai/ximi-y8y4).Given that at the end of June snow cover extension has not yet achieved its annual minimum size, it is seen how the cores of the glaciers are contained in boundaries identified as snow by MODIS.A similar view is shown in Fig. 4b.Here, we generated polygon-shaped glaciers from Terra and Aqua snow maps of 18 July 2012.Such areas are then overlapped to orthophotos taken about a month later.

Data
The Spectroradiometers on board Terra and Aqua satellites provide source data for the snow-covered area product (SCA) available in the National Snow and Ice Data Centre (NSIDC).Maps belong to Collection 5 (Riggs et al., 2006) of MODIS snow products, released after the reprocessing which began in September 2006, showing improvements in the cloud mask product and improved screening for erroneous snow.Cloud masks are derived from MOD35, with daily tempo-ral and 1 km spatial resolution.Each pixel not masked by clouds is processed by the snow-mapping algorithm Snomap (Hall et al., 2001).It consists of some tests and decision rules that identify snow in each pixel of a MODIS image, and it is able to detect snow even in dense forests (Klein et al., 1998).The procedure for snow detection employed by MODIS Snow Products works on a pixel-by-pixel basis involving a grouped-criteria technique which uses the Normalized Difference Snow Index (NDSI) and other spectral threshold tests (Hall et al., 1995(Hall et al., , 2001;;Klein et al., 1998).SCA is included in MOD10A1, a set of daily products available from February 2000 that contains data from the Terra satellite, while MYD10A1 are similar products derived from MODIS on board the Aqua satellite (launched in May 2002).The snow cover daily tile (SCA) contains a binary information (snow/not snow) for each pixel.MOD10A1 and MYD10A1 also provide daily snow albedo (Klein and Stroeve, 2002) and fractional snow cover (FSC) (Salomonson andAppel, 2004, 2006).Our choice of adopting the binary product instead of the fractional one is due to its wider use in the last decade, which has led to several assessments of its reliability (Hall and Riggs, 2007).MOD10A1 and MYD10A1 are generated through a mapping of pixels derived from the swath products MOD10L2 and MYD10L2 to their geographic locations in a MODIS-specific global sinusoidal projection.(Riggs et al., 2006) into three classes: snow, land (not snow) and cloud.Here, "Missing data" and "No decision" fields are treated equally to cloud pixels, where snow presence has to be investigated.We pre-processed daily Terra and Aqua SCA maps from January 2003 to December 2012.Except for a few particular days in which the recording failed, over this period two daily observations are available from Terra and Aqua.Accordingly, the cloud removal procedure requires as input two pre-processed snow cover maps per day and the DEM of the area.

Results
In this section we first provide some considerations about cloud obstruction and its negative effects on snow mapping in Alpine areas.Then, we show the results of the cloud removal methodology applied to the case study.Finally, the performance of the procedure is assessed in terms of its efficacy and accuracy.

Cloud obstruction over the Po River basin
The variable C d ranges from 0 to 1 over the year, depending on the weather.The histogram of its absolute frequency (0-366 days) is plotted in Fig. 5 for 2004.The distribution is based on the relationship between the spatial scale of atmospheric disturbances and the size of the basin.As can be seen, the most frequent values are those wherein the domain is covered by cloud for more than 90 % or for less than 20 %.
The mode of the distribution can be estimated around 95 %.
The peak is even more evident when analysing the Aqua images.The same happens for all the other years: it follows that in cloudy weather usually the Po Valley lies entirely beneath clouds, without information for snow mapping from MODIS.This is an interesting result since spatial filters such as step 3 estimate snow presence on the basis of data available in the same daily map.Thus, the distribution of C d suggests the introduction of some previous steps before the assessment of snow cover distribution over elevations.Processing MOD10A1 and MYD10A1 products of 10 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) we obtain that the mean annual cloudcovered part of the basin corresponds to a percentage between 44 and 56 % for the Terra MODIS and between 46 and 59 % for the Aqua MODIS.These values are slightly lower than those found by Parajka et al. (2006) for Austria.The higher percentages found in the Aqua images are due to the fact that Aqua crosses northern Italy a few hours later than Terra, in the afternoon, when cloud cover usually increases.In both cases these rates exceed 60 %, analysing only the Alpine regions located at elevations higher than 1000 or 2000 m a.s.l.(Table 3); such result is certainly related to the hours of the day when these satellites overpass the region (between the late morning and the afternoon).This confirms the evidence that in afternoon clouds usually persist, especially over mountain areas.The effect is partially related to cumulonimbuses, forming from water vapour carried by upward air currents, whose rising is hampered by orography.At the same time air circulation in the valleys may contribute in cleaning low altitudes from clouds.The influence of topography is particularly evident in Fig. 6: here, we plotted the number of cloudy days per pixel recorded by Terra MODIS in 2008 and 2011.It is possible to see that cloud pixels show patterns related to elevation (see Fig. 3a for the DEM).Similar patterns also occur in the other years, but are not reported here for brevity.
From Table 3 the values of average annual cloud cover fraction for different threshold altitudes confirm that the probability of being masked by clouds increases with altitude.Since uplands are also those where snow is expected in  winter and spring, the problem of cloud cover has a significant (and adverse) impact on the number of days in which the presence of snow can be detected.Average values of cloud cover can be calculated per quarter (or per month), in order to establish whether a seasonal trend exists (Table 4).While considering the entire domain JAS (July-August-September) is characterized by largely lower values of cloud obstruction, but the same trend does not apply for the part which lies above 1500 m a.s.l.In contrast, the mountain zone presents similar cloud cover percentages in the first, third and fourth quarters, but higher values are found in AMJ (April-May-June).

Performance of the cloud removal procedure
The performance of the procedure is assessed here.First, the efficacy of each step is investigated, which is its contribution in gap-filling over the cloud masks.Then, we focus on the accuracy of the steps and of their application in series.

Efficiency
Table 5 summarizes the mean annual percentage of cloud cover C a remaining after each step, proceeding from the  source Terra and Aqua products (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) to the final gap-filled map.The first step combines two "observed" daily maps (Terra and Aqua snow cover tiles): hence column "Step 1" represents the average annual fractions over which snow presence has been estimated using information from close days or cloud-free areas.This value is usually lower than 0.4, thanks to the well-known efficacy of combining Terra and Aqua data (Parajka and Blöschl, 2008;Gafurov and Bárdossy, 2009).
Step 2 is still based on recorded images, but maps are those of neighbouring days.This step, introduced by Gafurov and Bárdossy (2009), proves to be very effective, since after its intervention the remaining mean cloud cover generally decreases to less than 30 %. Step 3 involves the regional snow and land lines and results in an average decrease of the cloud mask of 3-5 %.This could seem a weak contribution; however, it should be considered that step 3 is involved only in days when clouds hide less than 50 % of the basin.Such condition is met on slightly more than half days per year; as a consequence, its contribution on these days is higher than it may seem looking at the average annual values.Moreover, the snow line is not used from June to September.From the monthly average, it is clear that the intervention of step 3 is not negligible in winter (Fig. 7).
The regional snow and land lines are shown in Fig. 8 for 2003 (Fig. 8a) and 2008 (Fig. 8b), during the period January-May.Also reported are the daily altitudes of the highest land pixel and of the lowest snow pixel, in order to appreciate the difference with a second approach for delineating snow and land lines (Gafurov and Bárdossy, 2009;Paudel and Andersen, 2011).It appears that the effectiveness of this step greatly increases when assessing such altitudes through average elevations: mean regional elevations of snow are every day lower than the altitude above which all pixels are snow.The same applies for land areas, whose average elevation is ordinarily higher than the lowest snow.In Fig. 8, the rising of snow average altitude is related to melting, while lowerings are caused by snowfalls over low elevations that bring down the mean values.Figure 8a shows the faster variability of the snow line in spring, when snow processes are fast and occasional snowfalls occur.On the other hand, the fluctuations of the maximum land elevation are always wider, especially in winter: this is due to the fact that such estimate can be highly influenced by anomalous local conditions, or misclassification errors from the source MODIS product.The trend of the minimum snow elevation shows that during most of the winter Terra and Aqua maps contain at least one snow cell around sea level.This could be result of specific local conditions occurring somewhere within the wide domain, as well as classification errors by MODIS Snomap (Hall et al., 2001).
Since in step 3 we consider four aspect classes while maintaining the regional scale of computation, we first investigated whether the "signal" of aspect exists even at the regional scale.Figure 9 shows the altitudes of the regional snow/land lines from January to the end of May (years 2005 and 2009).Lines are drawn per aspect classes (north: aspect > 315 • and aspect < 45 • ; south: 135 • < aspect < 225 • ).Daily average altitudes are computed only for days with less than 50 % cloud cover using output maps from step 2. It is seen that the impact of exposure exists even at the regional scale.Even though meteorologic variables may locally affect snow distribution on elevations, snow lines of the south class stand consistently above those of the north class.Except for days with extensive snowfalls within the basin, which lower the snow lines uniformly, snow line altitudes increase during the melting season and the south class maintains higher values for all the study period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012).Elevation differences can be quantified as of the order of 200-250 m as average value from January to May.As expected, such differences are less pronounced for the land line, which is mainly dominated by the great number of pixels at low altitudes covering the Po Valley.Overall, this demonstrates that a regional snow line approach which considers exposure improves the representation of topographic effects on snow distribution.Since aspect values can be derived directly from the DEM, our solution in step 3 introduces an improvement without any requirement of additional information.While the first three steps proceed in series (i.e. the output map of each is the input of the following), step 4 uses the most recent observation for each pixel (produced by step 1), as long as it falls within the previous week.Such operator is the busiest and reduces the mean annual values of cloud coverage to 2-4 %.Step 5 completes the cloud removal.Given the average annual cloud cover filtered by the last step (2-4 %), it seems to entail a little effect; yet its involvement is crucial, since it takes care of all those situations when cloud cover persists over the same area for several continuous days.On some dates the percentage of clouds remaining after step 4 are even higher than 10 %.
To highlight whether there are monthly differences in the effectiveness of each step, Table 6 contains   Step 2 leaves out a cloud coverage of 29 %.The intervention of step 3 affects 14 % of the domain and step 4 closes the gap-filling procedure.
cloud cover removed by the middle passages (e.g. the distance between a point of line "step 3" and a point of line "step 2" is the mean monthly contribution of step 3; the distance between a point of line "step 1" and a point of line "step 4" is the mean monthly contribution of steps 2 and 3).Line "step 4" is the mean percentage of cloud cover treated by step 5. Heavy seasonal patterns are difficult to spot; nevertheless two remarks seem appropriate.First, we note an expansion of the band between line "step 3" and line "step 2" during the winter.Such contribution represents the intervention of step 3; the enlargement stems from the fact that the regional snow line in winter is located at lower elevations and so enhances the number of cloud pixels lying over (then classified as snow).In contrast, the annual fluctuation of the regional land line is moderate, ruled by the wide lowland which experiences only occasional snow cover (Fig. 8).Accordingly, the transition elevations between the regional lines, where step 3 does not act, are restricted in winter.
A second consideration concerns the behaviour of step 5 (distance between line "step 4" and the x-axis): its average monthly contribution can exceed 10 % in periods characterized by very high cloud cover rates.The explanation lies in the fact that high monthly percentages normally entail several continuous days of widespread cloud cover, which reduce the probability of finding data in the previous week.A similar situation appeared for Austria in Parajka and Blöschl (2008), using 5-and 7-day temporal filters in the period 2003-2005.Thanks to this result, the importance of a step able to close the procedure is clearer.Figure 10 shows the result of the procedure run on 27 March 2003.Here, step 4 basically closes the procedure, thus the final product of step 5 is not shown.It is possible to appreciate changes in cloud distribution on the crossing times of Terra and Aqua: in this case the Aqua map does not provide additional data.Step 3 contributes in cleaning the very high and very low elevations, thanks to the regional lines, without affecting the transition band.The removed-cloud map of step 4 sheds light on snow distribution and its dependence on topography at the beginning of the spring.A second example is provided in Fig. 11 with reference to 22 January 2012.

Accuracy: 1 day additional cloud mask
The validation sample is composed of 25 days selected over the period 2003-2012.Results are shown in Table 7.The last two columns contain the performance of a basic 7-day backward temporal filter.The comparison is described in Sect.4.2.4.
Overall, our improved approach shows a good level of accuracy in all seasons.For most of the tests the degree of agreement D A reaches 96 %.A slight tendency to overestimate snow cover is evident; in spring, when melting rapidly reduces the snow-covered areas, the temporal filter in step 4 was expected to encompass this effect.However, overestimations occur also in autumn and winter.The maximum degree of disagreement (10 %) is obtained on 12 December 2012, after the application of a large cloud mask from 4 February 2012.Days before this date were characterized by a wide snowfall in the lowland, involving a later fast evolution of ground conditions.Considering the area above 200 m a.s.l., we get D A = 95 % for the same testing day.A clear underestimation affects the accuracy of the procedure on 9 October 2011.Looking at the source maps, this day appears as the first one with a well-developed snow cover on the entire mountain region, shrinking back in the next observations.The underestimation is mainly due to step 4, which finds a lack of snow in the previous period.In general, autumn is characterized by snow areas under development and any type of filter looking backward may underestimate snow.Tables 8  and 9 report the performance for each step of the procedure, relative to 3 March 2012 and 10 December 2004.The same analysis was carried out for the whole sample.Step 1 shows a high level of accuracy when outputs maps are compared to the original Terra products.We greatly limited its intervention by masking the source maps of both satellites: this measure was necessary to involve the other steps in the test.In general, more than 97 % of the cells classified by step 1 were correctly assigned to snow or land.Misclassifications may come either from changed ground conditions on the Aqua times of overpass or from the Aqua snow-mapping algorithm which uses MODIS band 7 instead of band 6 (Riggs et al., 2006).An example of poor performance for step 1 is shown in Table 9b for 10 December 2004.It was a cloudless day when snow cover had a wide distribution in the Alps.Terra provided a clear map with snow only on mountain reliefs.Aqua overpassed the region about 2 hours later and its binary product shows snow in the western plain.
Step 1 uses the latter to reclassify some additional clouds, contrasting with the source Terra product.As already found in Gafurov and Bárdossy (2009), step 2 has a high degree of accuracy.Even if it involves information from close days, the efficiency is ensured by the requirement of stationary ground conditions within the temporal window.We obtained a degree of agreement higher than 97 % for 21 days of validation, considering only cloud pixels processed by step 2. Table 10 shows the performance of the combined use of step 1 and step 2 both in terms of efficiency and accuracy.The performance of a simple 2-day backward temporal filter is reported for comparison.Each filter processes only a fraction of the additional cloud pixels N c,T , equal to A d,F /A d,T , where A d,F ≤ A d,T .Its degree of agreement is then evaluated on such processed cells.Values demonstrate that the conservative temporal filter (subscript ±2 days) is less effective while saving a definitely higher level of accuracy.As expected, the basic 2-day backward filter (subscript −2 d) provides a greater contribution to cloud removal.However, the absence of controls on what happens in days following a cloud observation decreases the accuracy, since in some cases the ground condition is changed.The choice of adopting a method based on a double control in step 2 ensures high-accuracy maps over which the following step can compute the regional snow line.
The regional snow and land lines defined in step 3 have shown a high efficiency whenever this method is involved in classifying a significant number of pixels.These elevations seem representative of the threshold altitudes above (resp.below) which all cloud pixels are likely snow (land).On av-erage, the percentage of cloud cover processed by this step is lower than that of step 2. Nonetheless the degree of agreement provided for the majority of the days is comparable.On dates when its operation involves more than 3 % of the additional clouds (10 cases), its degree of accuracy never drops below 96 %.An example is 11 November 2009, where 21 % of the additional mask is filled by step 3, and the individual value of D A reaches 99.2 %.On 8 December 2011 its contribution affects 12 % of the artificial clouds with an accuracy of 98.8 %.On 5 February 2004 it processes 14 % of the cover showing an individual D A equal to 99.4 %.Moreover we verified that the daily average altitude of snow areas is generally lower than the elevation of the highest land pixel, thus allowing a wider action of the regional snow line approach.
Performances of steps 4 and 5 should be examined in days when previous steps leave out a high percentage of cloud cover.In fact, the last steps are called to close the procedure solving all cases of doubt: critical situations are those where ground conditions change quickly, or within the transition band where the snow line method does not intervene.Therefore a lower accuracy was foreseen and must be accepted.In days when step 4 processes wide cloudy areas, the accuracy of the 6-day (backward) temporal filter is good.For example, on 31 March 2006 step 4 classifies by itself 89.5 % of the pixels artificially masked, providing a degree of agreement of 96.3 %.On 27 November 2004 its contribution involves 72.9 % of the additional cloudiness and classifies correctly the 94.6 % of them.A similar performance is provided on 17 May 2012 when its degree of agreement is about 97 % while removing 71 % of N c,T .
A relatively extensive operation by step 5 takes place on 10 December 2004 when it processes about 10 % of the additional mask.Here, the degree of agreement reaches 95 %.
Step 5 is involved even on 21 April 2004 and on 27 September 2006, processing about 5 % of N c,T with agreements of 95.0 and 99.8 %.This step has been validated even individually on a reduced sample composed of 6 days from 2004 and 4 days from 2009.In the setup for the individual assessment, it intervenes after step 1 removing the whole cloud mask.It provided an average degree of accuracy of 95.0 %, with a tendency to underrate snow cover in spring and to overestimate in autumn.A second setup of step 5 was tested individually.In this case, the first land observation after February and the first autumnal snow cover split the year in land and snow periods per pixel.This methodology, introduced by Gafurov and Bárdossy (2009), provides an overall degree of accuracy of 90.3 %.We are driven to think that the improvements we suggested for individuating such macro-seasons perform properly.In particular, basing the start flag for snow season on more snow observations, we protected the performance from wide overestimates in the late autumn.

Accuracy: spatio-temporal additional cloud mask
Through this supplementary validation we studied the performance of the procedure when applied to several subsequent days with significant cloud coverage.Four groups of clearsky days are chosen within the period November-May of years not yet involved in the validation (2003, 2007, 2008 and 2010).Cloud masks are derived from Terra and Aqua imagery extracted from the same seasons of other years.
Results are reported in Table 11.The average degree of agreement D A decreases from 95.7 to 94.4 %.The tendency to overestimate snow cover is still evident and affects the whole sample.In particular, a low degree of agreement is found between source and processed maps of 8, 9 and 10 February 2008.The misclassification is due to the 6-day backward filter of step 4. In fact, the previous week suffered a snowfall over mid and low altitudes.The filter uses such data for gap-filling after the artificial masks are applied, while in the source imagery snow cover disappears quickly in the lowlands.This outcome is confirmed even by Table 12, where we reported efficiency and accuracy following the series of step 1 and step 2. Here, A d,±2d /A d,T represents the fraction of artificial cloud mask filled by step 1 and step 2 while D A,±2d is the percentage of pixels correctly classified with respect to the clear-sky Terra map.In the first 4 days of validation for year 2008, step 2 is not very effective due to the invasive and continuous cloud cover introduced.Accordingly, step 3 is skipped given that cloud coverage does not reduce to less than 50 % of the domain after the first steps.Most of the clouds are reclassified by step 4.
Table 12 outlines the low effectiveness of the conservative temporal filter when cloud obstruction fills several continuous days.However, it confirms the good accuracy of the combined use of steps 1 and 2 in the processed areas.In the whole sample steps 1 and 2 reduce cloud coverage to less than 50 % of the domain only on 24, 25 and 26 December 2003.Thus, the employment of step 3 is seen at this point.The percentages of additional cloud pixels removed by step 3 are 20.8, 16.0 and 29.9 %, respectively.The related degrees of agreement are 96.0,98.2 and 99.1 %.Step 4 basically closes the procedure in all the days and it is responsible for most of the errors.In fact, the highest contribution of step 5 is observed on 26 December 2003 and affects 3 % of the additional cloud pixels.

Comparison with other cloud removal procedures
We assessed the improvements produced by this procedure in comparison with other existing approaches for cloud reduction.
First, to establish whether the stepwise procedure ensures any advantage compared to a mere backward temporal filter, we ran a 7-day filter on the same samples, applied to the merged Terra and Aqua images (outputs by step 1).This filter was tested in Parajka and Blöschl (2008) for Austria.Two sets of results were obtained.The first one is related to the percentages of cloud mask processed by each approach: while the new procedure treats the entire mask as in Gafurov and Bárdossy (2009), a filter limited to a preset time interval may leave out some cloud cover.Thus, the 7-day filter is expected to process a fraction of N c,T , equal to A d,7d /A d,T , where A d,7d ≤ A d,T .A second check concerns the reliability, evaluated as degree of accordance D A .Differences of accuracy between the new stepwise approach and the basic backward filter should reflect the improvements due to steps 3 and 5, since our procedure introduces temporal filters in steps 2 and 4.
Using the first validation strategy the comparison reveals that the stepwise procedure generates a concrete improvement when efficiency and accuracy are assessed in parallel (Table 7).Completing the classification, it processes even those pixels hidden by clouds for more than 7 days that are certainly the hardest.And yet, it does not entail a decrease in the accuracy.An actually higher performance of the 7-day filter emerges only on 23 May 2009, when its degree of accuracy D A,7d outweighs the other by 0.5 %.However, even in this case the 7-day filter leaves 1.3 % of the additional cloud mask.Conversely, on three validation days the accuracy of our procedure exceeds that of the temporal filter by more than 1.5 %.The percentages of additional clouds left out by the temporal filter are between 4 and 6 % for 5 days: focusing on our methodology, steps 2 and 4 require information on the previous days, as well as the 7-day filter.In such cases steps 3 and 5 come into play, allowing the classification of more pixels without entailing higher disagreement.Overall, the value of D A of the proposed procedure is 95.7 %.On the other hand, the temporal filter alone provides an average degree of agreement of 95.2 %.Moreover the variability in the accuracy of the performances is slightly lower in our procedure.
In the second assessment (Table 11) the improvement is lower due to the poor contributions of step 2 and 3.For most of the days the accuracy of the new method, D A , reduces to that of the basic 7-day backward filter.A slight gain emerges only on 24, 25 and 26 December 2003.Here, the procedure exceeds the performance of the temporal filter both in terms of effectiveness and accuracy by means of steps 3 and 5.
We also compared the procedure with respect to the one proposed by Gafurov and Bárdossy (2009).Steps 1 and 2 are equal in these two stepwise solutions.In Gafurov and Bárdossy (2009) the snow line and land line of step 3 are evaluated as the highest land pixel and the lowest snow pixel, if the study area is at least 70 % cloud free.The temporal filter is not involved, while two spatial filters compose steps 4 and 5, comparing each cloud cell with ground conditions of the nearest cloud-free pixels.The procedure is then completed with a "macro-cycle" approach working individually for each spatial unit, similar to our step 5 (for a complete description see Gafurov and Bárdossy, 2009).Differences lie in the way snow and land periods are evaluated.In Gafurov and Bárdossy (2009), the winter season coincides with the time window following the first snow observation, while land periods starts at the first land observation in spring.In contrast, our step 5 checks ground conditions on more days before identifying the snow-land cycles.Using the methodology by Gafurov and Bárdossy (2009), we obtained an overall D A of 91.8 % (σ = 7.2 %) and D A of 91.6 % (σ = 4.8 %) for the first and second validation strategy respectively.In most cases the two solutions provide similar results; however, the combination of steps 4 and 5 in our procedure ensures much higher performances in closing the removal.The identification of the "snow period" through the first observation, in autumn or winter, produces excessive overestimates for this case study.Micro-cycles are very common before winter, and the fluctuations of the area covered by snow can be very wide in a short time window.Moreover we found a very low effectiveness of spatial filters which work at the pixel scale, giving slight contributions over all the seasons.

Discussion and conclusions
In this work we investigated the use of MODIS products for snow mapping over the Po River basin.The large-scale analysis, required when one deals with the major European rivers, often suggests the use of imagery from remote sensors.Daily snow cover maps for the whole domain require the adoption of a cloud removal procedure able to estimate ground conditions beneath clouds, in order to preserve the temporal resolution of the data.
The first question we tried to answer is: how much does cloud obstruction affect MODIS daily snow products?Looking at our study region, this influence seems certainly more disturbing than what we expected.The average annual values of cloud cover range from 46 to 59 % of the basin for Aqua observations and from 44 to 56 % for Terra.These percentages are slightly lower than those found in Austria (Parajka and Blöschl, 2008).However, seasonal trends show that in the Alpine space this disturbance increases during the melting season.Focusing on spring, for most of the years (from 2003 to 2012) average values of cloud cover overcome the 70 % of the region which lies above 1000 m a.s.l.Spring is the most interesting period for hydrologic studies since snow dynamics are faster and daily observations allow us to follow gradually the patterns of melting, determining the contribution to streamflows.For this season, our analysis brings into question the real benefit of maps with temporal resolution of the order of 1 day.Causes must be sought in the times of the two satellite overpasses, that fall between the late morning and the afternoon.Further investigation may be directed to a thorough identification of factors that determine seasonal and spatial patterns of cloud obstruction in MODIS snow maps.
Results obtained from the analysis of cloud obstruction justify the use of a cloud removal procedure.Several approaches have been proposed to reduce negative impacts of clouds on snow mapping with high temporal resolution.After some first attempts we noticed that methods already applied to other regions introduce assumptions that did not fit well for the climatic and topographic features of the Alps (Gafurov and Bárdossy, 2009;Gao et al., 2010;Paudel and Andersen, 2011).Northern Italy is mainly characterized by a range of elevations where snow undergoes several short cycles of accumulation and melting, thereby complicating any estimation of what happens under clouds.Thus, we tried to combine and improve steps from previous studies, in order to maximize the performance of the procedure in the considered case study.Our procedure requires only MODIS binary snow maps and a DEM of the domain, without additional data from other satellites nor information on land uses.It entails a computational effort slightly higher than a simple temporal filter, while completing the gap-filling even when clouds hide the ground for extended time windows.We proposed two validation strategies (1 day and multi-day) for evaluating the reliability of the methods.Both are based on the introduction of fictitious clouds on clear-sky images before running the code (Gafurov and Bárdossy, 2009).The cloud removal procedure ensured an average degree of accuracy of 95.7 and 94.4 % respectively.A slight tendency to overestimate snow areas was detected.The performance of the procedure was compared to that of a basic 7-day backward temporal filter in terms of efficiency and accuracy.Improvements in the accuracy are mainly due to the regional snow/land line step (Parajka et al., 2010).In this step we introduced the calculation of snow and land line altitudes per aspect classes, demonstrating that the signal of exposure exists even at the regional scale.When snow covers are well developed (mainly in winter and spring), such method plays a leading role.The efficacy of the novel procedure increases by means of step 5.In fact, it involves an improved yearly cycle seasonal filter, which can classified pixels regardless of the amount of data existing in the nearest days.It is based on the definition of two periods per year and per pixel, one affected by snow cover and the second without snow on the ground.The accuracy of this final step proved to be only slightly lower than the others.
We also want to highlight the good results provided by a simple temporal filter limited to 7 days, which uses the latest information available for each cell (overall accuracy of 95.2 and 94.2 %).Here, the greatest benefit is the ease of use.The limitation is that any restricted temporal filter leaves out a percentage of cloud cover which can be significant in the Alps, where clouds may cover the same area for several continuous days.
Approaches directed to the individuation of a typical condition for each cell on the basis of the period of the year (Gafurov and Bárdossy, 2009;Paudel and Andersen, 2011) seem less suited to this environment.The identification of the first snow observation as the starting time of a window when a certain location is expected to be snow covered leads to excessive overestimates.We verified that short-term snow covers are common over the transition altitudes, both in autumn and spring.
The final cloud-free maps provided by our procedure appear to be a reliable starting point from which it will be possible to derive information on snow resource for the Po Basin.The integration with snow depth measurements and snow density estimations will provide assessments on SWE amounts and fluctuations over the last 10 years.Moreover, a daily binary information on snow presence is in itself a useful tool for many applications at the regional scale.For example, snow duration and snow extent can be directly derived and associated with the inter-annual variations of Alpine temperatures and precipitation for studying snow cover climatology.
The Supplement related to this article is available online at doi:10.5194/hess-11-4579-2014-supplement.

Figure 1 .
Figure 1.Flowchart of the cloud removal procedure representing the step order and the temporal window considered in each stage.

Table 2 .
Thresholds of the number of snow (land) observations (following the first one) required for identifying the beginning of the snow (land) season in step 5.The threshold is set per altitude ranges.Altitude range n s n l [m a.s.l.]

Figure 2 .
Figure 2. Criterion for identifying the start of the land season.The same rule is applied to the snow-start flag, in autumn or winter.More land/snow observations are needed in spring/autumn to split the year into two periods.The black line represents the series of observations for a reference pixel across the melting season.

Figure 3 .
Figure 3. (a) DEM of Italy at 500 m spatial resolution with the contour line of the Po River basin.(b) Hypsographic curve of the Po River basin derived from the DEM.
3.2), clearsky days provide an indicative value of the fluctuation of the snow-covered area throughout the year.For instance on 1 January 2005 snow is about 30 % of the clear-sky pixels, while this percentage decreases to 5 % at the beginning of May (average elevation of the snow-covered pixels of 2640 m a.s.l.) and less than 0.5 % at the end of June, when the mean altitude of the snow-covered cells exceeds 3260 m a.s.l.Analysing the most clear days without snowfalls in the lowlands of the 2004/2005 and 2005/2006 winter months, values between 20 and 30 % seem a representative range for the long-lasting snow-covered fraction of the basin.Notice that the area-elevation curve of the Po Valley shows that about 30 % of the area lies above 1000 m a.s.l. and 20 % above 1500 m a.s.l.To emphasize

Figure 4 .
Figure 4. (a) Comparison between the digital map of Lombardy's glaciers (sketched in blue) and the snow distribution (red) from MYD10A1 at the end of June 2005 (source Regione Lombardia https://www.dati.lombardia.it/Territorio/Ghiacciai-e-Nevai/ximi-y8y4).(b) Snow polygon derived from MOD10A1 of 18 July 2012 overlapped on the orthophoto AGEA 2012 of Lombardia.The latter is a composite image from air flights operated at the beginning of September 2012 (source Regione Lombardia http://www.cartografia.regione.lombadia.it).

Figure 5 .
Figure 5. Absolute frequency of the daily cloud percentage C d from Terra MOD10A1: year 2004.

Figure 7 .
Figure 7. C m , average monthly values of C d remaining after the implementation of each step for year 2007 (a) and 2009 (b).

Figure 8 .
Figure 8. Regional snow and land lines from January to May in comparison with the elevations of the highest land pixel and the lowest snow pixel, year 2003 (a) and 2008 (b).Values are provided only when output maps by step 2 are less than 50 % cloud covered.
average cloud cover fractions per month C m remaining step by step.The values refer to 2010. Figure 7 shows the same results plotted for 2007 (a) and 2009 (b).The vertical distance between points of different steps represents the average values of

Figure 10 .
Figure 10.Cloud removal procedure applied to 27 March 2003.Aqua and Terra maps have a cloud coverage of 98 and 83 %, respectively.Step 2 leaves out a cloud coverage of 29 %.The intervention of step 3 affects 14 % of the domain and step 4 closes the gap-filling procedure.

Figure 11 .
Figure 11.Cloud removal procedure applied to 22 January 2012.Aqua and Terra maps have a cloud coverage of 77 and 91 %, respectively.The cloud removal is basically completed by step 4.

Table 1 .
Intervention rules for step 2. Cloud pixels to be reclassified are those at day d.The table displays only those combinations of observations between d − 2 and d + 2 that admit the intervention of step 2. In all the other possible cases this conservative filter is skipped.

Table 3 .
C a (Z > z T ): mean annual cloud cover fraction for different altitude thresholds over the Po River basin.Daily snow products MOD10A1 and MYD10A1 are downloaded from NASA's EOSDIS "Reverb", a metadata and service discovery tool (web portal: http://reverb.echo.nasa.gov).Snow cover daily tiles are provided in sinusoidal projection.Tiles are 10 • by 10 • at the equator.The study area is entirely covered by the h18v04.Snow maps of Northern Italy are then projected into the UTM WGS84 coordinate system, Zone 32 band T, and overlapped with a mask of the Po River basin.We reclassified MODIS snow cover maps from the original pixel classes z T = 0 m a.s.l.z T = 1000 m a.s.l.z T = 2000 m a.s.l.

Table 4 .
Cloud cover fraction over the Po River basin averaged per quarter from MOD10A1.Values refer to the entire basin and to the area higher than 1500 m a.s.l.

Table 5 .
Mean annual cloud cover fraction remaining after each step for the Po River basin.

Table 6 .
Average monthly cloud cover fraction remaining after each step for year 2010.

Table 7 .
Overall accuracy of the cloud removal methodology tested on several days and comparison with a basic 7-day backward temporal filter.Values in [%].

Table 8 .
Cloud, land and snow percentages (C d , L d , S d ) in clearsky maps of 3 March 2012 and 10 December 2004, and percentages after the application of the additional cloud cover (A d ).

Table 9 .
Results of the validation step by step for 3 March 2012 and 10 December 2004.

Table 10 .
Efficiency and accuracy of the conservative filter of step 2 (subscript ±2 d) compared to a basic 2 day backward temporal filter (subscript −2 d).Both filters are applied to maps preprocessed by step 1.Values in [%].A d,T A d,-2d /A d,T D A,-2d A d,+/-2d /A d,T D A,+/-2d

Table 11 .
Overall accuracy of the cloud removal methodology tested using spatio-temporal additional cloud masks (MOD10A1 SCA of 24 December 2003 is replaced with MYD10A1 SCA since data are not available).Values are compared to the performances of a basic 7-day temporal filter.Values in [%].

Table 12 .
Efficiency and accuracy of the combined use of steps 1 and 2. The test is performed using spatio-temporal additional cloud masks.