Comparison of MODIS and Model-Derived Snow-Covered Areas: Impact of Land Use and Solar Illumination Conditions

: Moderate resolution imaging spectroradiometry (MODIS) snow cover accuracy has been assessed in the past at different scales, with various approaches and in relation to the many factors influencing the remote observation of snow-covered areas (SCA). However, the challenge of fully characterizing MODIS accuracy over forest sites is still open. In this study, we exploit 5 years of data from the upper river Adige basin at Ponte Adige (Eastern Italian Alps) to condition an enhanced temperature index snowpack model accounting for model parameter uncertainty by using the Generalized Likelihood Uncertainty Estimation (GLUE) methodology. The simulated SCA is then compared with MODIS retrievals through a range of different statistical metrics to investigate how land use and solar illumination conditions affect such comparison. In particular, the Overall Accuracy index (OA) is used to quantify the agreement between satellite-derived and simulated SCA on a pixel-by-pixel basis. Analyzing the spatial variability either of the median OA and its range shows that illumination conditions over forested canopies represent a major source of uncertainty in MODIS SCA. Exploiting this finding, we identify the minimum level of incoming short-wave radiation for accurate use of MODIS SCA in forest areas.


Introduction
Understanding and predicting snowpack dynamics is of crucial importance to characterize the water cycle in mountainous regions [1]. Snowpack prediction models are an essential tool for water resources management activities, such as hydropower energy production planning, irrigation and providing early flood warnings. Approaches to snowpack prediction range from empirical models (e.g., simple temperature index models) to more sophisticated physically based energy-balance models [2,3]. The advantages of an energy balance model are that it, in theory, requires less model calibration and has the potential to provide forecasts that account for climate or land cover change [4]. However, using an energy balance model may require more data, and model performance relies on the availability and quality of the additional needed climate input data. Prior studies comparing temperature index and energy balance models have not been conclusive as to whether one approach is better than the other [5][6][7].
Over the last two decades, a class of models which exploit both air temperature and solar radiation in the snowmelt computation, termed enhanced temperature index models (ETI), has shown to represent a viable approach for snowpack predictions [8]. ETI models are recognized to provide a promising approach to modeling the snowpack at the catchment scale with fewer input data than energy-balance models, but allowing better model parameter transferability than standard temperature index [8,9].
Remote sensing of snow cover is increasingly exploited for data assimilation [10,11] and assessment of snowpack models' performance [12,13]. Among the different remote sensing products, snow-covered area (SCA) maps retrieved from moderate resolution imaging spectroradiometry (MODIS) are largely used, due to their frequent temporal availability (1 day) and their relatively high spatial resolution (500 m). However, like all optical-based satellite sensors, accuracy of MODIS products is influenced by atmospheric conditions like overcast sky [14][15][16][17][18], land use, especially whether forested areas are present [14,15], and snowpack conditions such as shallow or patchy snow [19].
Such issues in MODIS snow detection accuracy might hamper the effectiveness of product integration within different model frameworks. In particular, remote sensing of snow in forested sites is challenging, as pointed out by several works carried out comparing MODIS-based SCA with simulations of snowpack models either based on the energy balance approach [17] or using an ETI model [20]. Outcomes of such studies show a clear drop in the average agreement between modeled SCA estimates and MODIS retrievals related to the forest sites distribution. In other words, regardless of the modeling framework, the forest canopy obstructs a reliable validation of simulated SCA using satellite-derived products.
Indeed, detecting snow under forest canopy by means of remote sensing techniques presents several issues related to the obscuration introduced by the forest canopy, which makes extremely challenging a reliable remote sensed detection of the snow accumulated underneath the trees. In cases of moderate canopy cover, up to around 60%, optical sensors can identify snow in the gaps between the trees [21]. Thicker canopies can hide the snow entirely, due to the heavy shadows. This error source depends on canopy structure, topography and viewing geometry. In particular, the forest effect is exacerbated by wide-view angles of MODIS as the viewable gap fractions through forests decrease when the view zenith angles increase [22]. Not surprisingly, contrasting results are reported in the literature. Most of the studies show that MODIS generally underestimates SCA in forested sites [14,15,21]. However, results reported by [14], for a small mountainous catchment in Northern Slovakia, show that the combination of Terra and Aqua MODIS images can enhance the overall accuracy of SCA prediction compared with extensive snow measurements. This procedure leads to good performances even in forested areas, with values of overall accuracy greater than the ones found for other land covers. The latter one partially agrees with other studies showing relatively good performances of MODIS even on forest sites, even though is still reported that forest represents a major issue affecting MODIS retrievals' reliability [23,24].
However, the abovementioned studies did not go through how MODIS snow detection behaves within forest sites, but they focus on MODIS reliability in dividing the study area into forest sites and open sites, mainly for investigating the land-use impact on MODIS retrievals. The research effort carried out in this paper aims to fill this gap, explicitly investigating MODIS snow detection performances on forested areas, comparing MODIS-based SCA with predictions obtained from a recently proposed topography-based distribution function snowpack model, TOPMELT [20].
The following research questions are addressed: (i) analysis of seasonal accuracy of MODIS retrievals on forest sites and comparison with other land uses; (ii) designing of a novel methodology to discriminate forest areas over which the accuracy of MODIS data is still deemed acceptable and therefore is worth using for model validation or data assimilation purposes.
The paper is organized as follows. Section 2 provides information about the case study area as well as TOPMELT theoretical basis and uncertainty analysis. Then, in Section 3, the results from the comparison between satellite-derived and simulated SCA are shown, and they are discussed in Section 4. Finally, Section 5 reports some concluding remarks as well as suggestions for future developments.

Study area and dataset
The study area considered in this work is the Adige river basin at Ponte Adige, in the Eastern Italian Alps (Figure 1). The basin area is 2179 km 2 with an average elevation of about 1900 m, ranging from 230 to 3860 m.
Precipitation in the study area is controlled by two main gradients: i) a West-East general increase of mean annual precipitation, ranging from 490 mm/yr at the valley bottom in the western and drier portion of the basin to 800 mm/yr in the eastern portion; ii) an orographic enhancement of the precipitation, with a lapse rate which ranges from 122 to 445 mm/km.
Typically, at higher elevations (over 2000 m), early snowfalls occur during October, and by the beginning of November, a permanent snowpack can be observed, which usually lasts until August. At intermediate elevations (from 1000 to 2000 m), snowfalls mostly occur between December and March, while for the rest of the year, the precipitation can be classified as rainfall. Figure 1b shows the different land uses inside the catchment. Land cover features were defined based on the well-known Corine dataset by aggregation into 6 classes (urban, agriculture, forest, grassland, bare soil and water bodies). Table 1 provides a summary of the aggregated land classes in terms of percentage area and mean elevation. Moreover, about 2% of the catchment area is covered by glacier. Glaciers were removed from the analysis because they might be misclassified as snow by MODIS sensor [23]. Forested area is mainly evergreen and located at low-intermediate altitudes.
Half of this area ranges between about 1200 m and 1850 m with an averaged value of 1514 m, while only about 10% has an elevation greater than 2000 m. Grassland dominates intermediate-high elevations, and the bare soil and rocky up-crops class covers the highest part of the catchment. Water bodies and urban areas cover a negligible part of the basin, i.e., 1.75% of the total area. Meteorological data are available from 42 weather stations installed over the case study catchment, all equipped with rain gage, recording hourly precipitation and temperature data. Besides, daily measurements of snow depth are available at 28 snow stations spread across the area (Figure 1c). The dataset also includes snow-covered-area maps, retrieved from MODIS sensor imagery. These maps are available daily at a spatial resolution of 250 m [4,5]

MODIS Snow Detection Algorithm
The MODIS snow cover maps exploited in this work (referred to as MODIS-EURAC) are based on the methodology described in [16,23] specifically developed for use in mountain areas. This algorithm allows the generation of snow maps with improved resolution of 250-m resolution with respect to the standard MODIS standard product (MOD10A1, MYD10A1), which have 500-m ground resolution. The algorithm was extensively validated by comparison with high-resolution snow maps derived from Landsat 7 ETM+ images and with snow depth measured from ground stations in selected test sites in Austria, Slovakia, Germany and Italy [16]. Moreover, a comparison with the NASA standard SCA products MOD10 (MYD10) was performed. The comparison of the proposed algorithm with both Landsat and MOD10 (MYD10) was conducted over the same overlapped area, which indicated an accuracy in the ranges of 88.1%-93.6% and 85.4%-95.3%, respectively. The lower accuracy values were found in forested sites, while in open free areas, the performances of the three products were quite similar, but MODIS-EURAC has the advantage of detecting more detailed features compared to the 500 m MOD10 (MYD10) maps. For the comparison with in situ measurements, the detected accuracy ranges from 94% to around 82%, with the lowest accuracies found in very rugged terrain. Moreover, the accuracy is lower for stations located in north-facing slopes, which may have a limited sun illumination in winter during the early morning acquisition. Another source of error is the misclassification with clouds, especially during summer [16,25] when patchy, dirty snow cover occurs together with clouds formation because of convective storms.
The major drawback of optical sensors is the lack of penetration through clouds, which hampers the acquisition of data during overcast days [10]. To address this issue, a simplified approach is adopted, leaving out from the data sample the days where MODIS detected a cloud coverage greater than 10%, following [26]. This criterion applied to the considered period (01 Jan. 2011 to 30 Sept. 2016) selected 550 maps, deemed to be enough to perform a consistent comparison between TOPMELT and MODIS-based SCA.

Snowpack Model: TOPMELT
Snowpack is simulated by using TOPMELT [20], an enhanced, semi-distributed temperatureindex model which uses air temperature and solar radiation to compute snowmelt. TOPMELT exploits a statistical representation of the distribution of clear-sky potential solar radiation by discretizing the full spatial distribution of solar radiation into a number of radiation classes for each elevation band. Thus, the computation required to generate a spatially distributed water equivalent reduces to a single calculation for each radiation class. This turns into a potentially significant advantage when parameter sensitivity and uncertainty estimation procedures are carried out. A detailed description of the model is given in [18]; thus, only a summary is reported here.
In TOPMELT, the basin area is subdivided into elevation bands to account for air temperature variability with elevation. Then, each elevation band is subdivided into a number of radiation classes. This is carried out by dividing each elevation band into a number of equally distributed radiation classes, where the i-th class contains the band sub-area corresponding to the i-th percentile of the incoming solar radiation energy. Therefore, the model spatial domain is represented by nb elevation bands and by nc radiation classes for each elevation band. Clear-sky short-wave solar radiation is computed at each element of the Digital Terrain Model (DTM) by taking into account topographic shadowing and complex topography. Since the model uses radiation values averaged over a given time interval, maps of potential radiation averaged over time are also computed. The spatial distribution of time-averaged clear-sky solar radiation is calculated over each elevation band, and equally distributed radiation classes are identified. For each radiation class, the mean daily cumulated clear-sky radiation value is computed (termed Radiation Index RI herein; (MJ m −2 d −1 )) and used in the snowmelt computation. Radiation is pre-processed and provided as model input.
Snow accumulation is computed based on air temperature and precipitation data from the available weather stations. Mean areal precipitation estimation is based on use of Thiessen's polygons, whereas air temperature data are used to estimate a spatially uniform dynamic vertical lapse rate over the whole basin. Every hour of available temperature records is linearly interpolated in order to estimate the lapse rate, which varies hour-by-hour accordingly.
Such mean areal precipitation can be corrected by means of a Precipitation Correction Factor (PCF) in order to account for the possible poor spatial representativeness of the available snow stations. Then, TOPMELT redistributes the sub-basin averaged precipitation, estimating the precipitation occurring at each elevation band, Pi (mm h -1 ) and introducing a vertical gradient G (% mm km −1 ) in order to take into account the orographic convection effect that leads to a larger amount of precipitation at higher altitudes. The total precipitation Pi is then classified as snowfall if the elevation band temperature Ti is lower than a threshold temperature Tc.
For the generic model cell given by the i-th elevation band and the j-th radiation class, day hours snowmelt rate , ( ) (mm h −1 ) at time t is computed, taking into account air temperature, clear-sky solar radiation and albedo, as follows: where: ( ) is the temperature of the i-th elevation band; is the class radiation index; [mm m 2 °C −1 W −1 ] is the combined melt factor, accounting for both thermal and radiative effects; ( ) [-] is the albedo of snow, = 0 °C is a threshold base temperature. Snow albedo is assumed to decrease over time starting from an initial value typical of fresh snow, albs. Following [27], the snow albedo is estimated through this relation: where β2 is a dimensionless parameter and the sum of Ti represents the sum of the positive hourly temperatures exceeding the threshold temperature Tb since the last snowfall (i.e., k hours). During night hours, snowmelt is estimated by accounting for air temperature only, while during rain-onsnow events, melting is computed depending on air temperature and on the energy provided by rain [28]. For each model class, the snow water equivalent ( , ) is updated by accounting for snow accumulation, rain-on-snow, melt and freezing water. Water due to snowmelt or rainfall is first retained in the snowpack as interstitial water ( , ) that propagates through the snowpack when exceeding the threshold of 10% of water equivalent. When air temperature is less than the threshold base temperature, part of the liquid water refreezes, and it is added to the current water equivalent.
It is worth noting that the described approach, based on aggregating pixels with similar solar [20] radiation, allows for high computational efficiency, typical of the simpler temperature index model ensuring at the same time the degree of accuracy and the more physically based approach of ETI models. Accordingly, TOPMELT might be an appealing tool for Montecarlo-like analysis, such as uncertainty estimation.
In principle, different setups of TOPMELT can be built, depending on the number of elevation bands nb (i.e., the width of each band), the radiation classes nc and the time span over which solar radiation is aggregated. Following [20], in this work, the width of the elevation bands was set equal to 100 m, each of them divided into 10 classes of solar radiation. The latter one was aggregated over a period of 4 weeks. A sensitivity analysis carried out by [20] showed that such TOPMELT setup applied on an Alpine catchment in South Tyrol guarantees good accuracy in terms of SCA simulations.
In addition, TOPMELT includes a vegetation module, which considers the transmission and attenuation of radiation through a forest canopy, precipitation interception and unloading, snowmelt and sublimation of intercepted snow.

Snow Compaction Model
In addition, to calculate the actual snow depth from the computed water equivalent, TOPMELT was integrated with a snow compaction model, following the scheme by [29], which computes snow density evolution over time according to a viscous compaction mechanism, driven by the combined effect of gravity and temperature. The relative density of freshly fallen snow accumulating over the existing snowpack is estimated for the i-th elevation band according to [30]. The relative compaction is calculated over time for the i-th elevation band and the j-th radiation class, treating snow as a viscous media compressing under the average weight of the bulk layer. The snow viscosity is calculated accounting for the effects of the liquid water content, snow density and temperature, exploiting the simplified model of [29] proposed by [30].

Spatial Representation of the Snowpack
Because of its distributed approach, TOPMELT can simulate the evolution in time of the spatial distribution of the SWE. Therefore, it can be validated against both distributed data (e.g., MODIS maps) and snow depth observations in order to check simulated snow-covered area and punctual SWE simulations, respectively. The SCA is calculated considering a pixel as snow-covered when the simulated SWE exceeds a threshold of 5 mm SWE, adopted by [31] to compare MODIS maps and ground snow observations across Europe. Moreover, a preliminary analysis was performed estimating the fractional snow cover area over a 2-year period (01 Oct 2012 to 30 Sept 2014) by means of TOPMELT using a set of different threshold (i.e., 0, 5, 10 and 15 mm SWE) to consider whether a pixel is snow-covered or not. However, no particular differences between the fractional snow cover area pattern were found, regardless of the threshold adopted, according to other similar studies [17].  Figure 3 shows an example of a categorical pixel-by-pixel comparison between MODIS-based and TOPMELT SCA relative to 25 Apr 2013. The categorical metrics derived from the contingency table were computed (i.e., "hits", "false alarm", "correct rejections" and "misses"), resulting in a generally very good agreement between what MODIS detected and what TOPMELT simulated. However, few "misses" (MODIS detects snow while TOPMELT does not) and "false alarm" (MODIS does not detect snow while TOPMELT does) are present, mostly in between "hits" areas (both MODIS and TOPMELT detect snow) and "correct rejections" areas (both MODIS and TOPMELT do not detect snow).

Model Parameter Conditioning
TOPMELT model parameters were conditioned applying the Generalized Likelihood Uncertainty Estimation (GLUE) procedure [32], comparing the simulated snow height with the records of 28 snow stations spread across the catchment (see Figure 2), over the period 01 Oct 2012 to 30 Sept 2014. During this period the precipitation regime shows significant snowfalls and therefore periods of high accumulation and melting dynamic. It is worth noting that summer months (July, August and September) were left out from this period since no snow was observed across most of the ground stations. A sample consisting of a 1000 parameter set was generated sampling the parameter space through the Latin hypercube scheme [33]. A preliminary sensitivity analysis of the model permitted to select the most sensitive parameters subject to random sampling ( Table 2) as well as their ranges. Table 2 summarizes the selected parameters as well as the respective range of variation. The value of the other TOPMELT parameters, not reported in Table 2, were set accordingly to a parameter optimization technique. In particular, the Snow Correction Factor (SCF) is set equal to 1.4, in agreement with those reported in previous works for the Alps [34].
Acceptable parameter set (i.e., behaviorals) selection requires the definition of a likelihood function and a threshold value for such a likelihood to discriminate between behavioral and nonbehavioral simulations. Selection of both likelihood function and thresholds are subjective choices of the modeler [35][36][37]. In this case, the likelihood function is the Nash-Sutcliffe efficiency (NS) index (see Equation (3)) evaluated at each of the snow stations available, as follows: where SDsim is the simulated snow depth at time step i, SDobs is the observed snow depth and SDobs,mean is the averaged observed snow depth over the simulated period. N is the total number of time steps and is equal to 108. For each parameter set, a set of NS values equal to the snow station number (amounting to 28) is obtained, and the median of the NS values is computed. A specific parameter set is deemed acceptable if it leads to a median NS exceeding a threshold set equal to 0.9. Therefore, behavioral identification is based on a multi-site approach, adopted to account for the significant spatial variability affecting snowpack dynamics [18]. The retained parameter set (630 out of a population of 1000) is used to simulate the predictive uncertainty range of SCA and perform the comparison with MODIS-based SCA maps. For this, the simulated fractional snow cover area was compared with the MODIS-based one in terms of Nash-Sutcliffe efficiency index over the period ranging from 01 October 2012 to 30 September 2014. Figure 4 shows the median NS index computed at each station, estimated comparing snow height observations with the median of the simulated snow depth carried out tuning the parameters shown in Table 2. Across the stations, the performance is generally good, with a median NS index equal to 0.91, consistent with what was estimated by similar comparisons carried out in other works [3,38]. An NS index ranging from 0.90 to 0.99 applying the enhanced temperature index model developed by [28] was found on a US catchment by [3], while [38] used the SWAT snow module in the Upper Adige river basin, calibrating its snow-related parameters with a two-year calibration period, resulting in an average NS index of 0.80.

Snowpack Model Validation
Furthermore, visual inspection of simulated and observed snow height showed that the snow duration is also well-reproduced by the model. Both median NS values are significantly lower than the median NS estimated across the stations. This happens since simulated snow height overestimates the observed data, especially from January to March 2014 where extremely high snowfall occurred. Pavicolo station is placed at 1400 m, while Meltina is at 1130 m, both having the southernmost aspect among all stations considered, about 260 degrees, counted counterclockwise from East. This southern aspect might have caused higher local temperatures compared with other locations at the same elevation, which were not well reproduced by the simplified approach of TOPMELT. Indeed, TOPMELT assumes uniform temperature at a given altitude, and this might have led to an overestimation of the solid precipitation that explains the bias between simulated and observed snow height during months.
Despite the accuracy drop that occurred at the abovementioned stations, no significant pattern was observed related to topography or elevation affecting the comparison between observed and simulated snow height across the catchment.
Finally, a seasonal estimation of NS index was performed, resulting in a median NS index during snow accumulation months (from December to March) equal to 0.89, while in spring and summer it rises up to 0.99, since most of the snow stations are snow-free. Again, no significant trends related to either the location or the height were found in winter nor in summer.
This means that TOPMELT is a suitable tool to simulate snowpack distribution, in particular the snow cover extent, and its related complex spatial variability over the catchment of interest.

Satellite-derived and Simulated Snow Cover Area Comparison
The TOPMELT-based fractional Snow-Covered Area (fSCA) and its relative parameter uncertainty estimated using the behavioral 630 parameter set is compared to the MODIS observations for the simulated period ( Figure 5). fSCA is defined as the ratio between the number of pixels classified as snow-covered and the total number of pixels excluding clouds and glaciers, as shown in Equation (4).
where Nsnow is the number of snow cover pixels according to MODIS dataset or TOPMELT estimates, Ntot is the total number of pixels representing the overall catchment area and Nclouds are the pixels classified as cloud by MODIS. Fractional SCA (fSCA) time series shown in Figure 5a are estimated considering the overall catchment area, i.e., regardless of different land uses. fSCA during both winters (December to March) included in the considered period is overestimated by TOPMELT, while the first part of melting season (April to June) is reasonably well reproduced. During summer months (July to September), MODIS fSCA ranges between 10% and 20%, while TOPMELT simulates a basically null fSCA. In October 2012, after a snowfall event, some days with temperatures above the seasonal average were observed. While the agreement between observation and simulation is good in terms of fSCA, immediately after the snowfall (fSCA about 65%), the melting phase is largely amplified by the model, leading to an underestimation of fSCA, likely because of such temperature anomalies. Figure 5b shows the comparison between fractional SCA considering only forest-free areas in order to analyze the effect of the canopy on the results, following [15]. This means that, based on the estimates provided in Table 1, a fraction of 32.17% of the catchment is neglected. A substantial increase of the agreement between MODIS-based and TOPMELT fSCA is observed during winter months for forest-free areas compared to the whole catchment case, while some disagreement arises during summer months. An in-depth discussion on this point is reported later in this section.
On the contrary, Figure 5c reports the same comparison but considering only forested areas. According to the previous findings, a clear disagreement stands out during winter months, when MODIS-based fSCA is always significantly smaller than the TOPMELT one. In summer months, the agreement is almost perfect, with vanishing fSCA for both sources during this period. Table 3 reports the season-averaged fSCA mean absolute error (MAE), estimated respectively on the overall catchment, the forest-free area and the forested sites. During the accumulation months, the mean absolute error over forested areas (MAE = 30.2%) is higher than the other two cases (Table 3), supporting the general findings reported in the literature concerning the effect introduced by the forest canopy. Moreover, the period 2012−2014 (i.e., the period over which MAE was estimated) was the wettest of the available time window (2011-2016) as shown in Figure 2. Since satellite-derived products tend to underestimate fSCA under forest canopy, the snowfall amount above the average that occurred during 2012−2014 might have further enhanced this feedback, leading to such high MAE. A similar comparison was carried out on a Norwegian catchment [35], resulting in an fSCA RMSE ranging between 6% and 22%, even though the authors considered open and forest sites together, without estimated the RMSE for each land cover.
The abovementioned findings support the idea that the presence of trees weakens the spectral signature from the underneath surface [17], hampering an optimal detection of snow cover by remote sensing devices. Among others, [39] argued the same conclusion comparing MODIS-based fractional SCA with the simulations of an energy budget snow model (SPH-AV). Likewise, [18] conducted an extensive MODIS validation in Austria, reporting low MODIS accuracy, evaluated by means of insitu snow depth observations, from November to March in forest areas. We should also take into account that most of the forest in the study catchment is evergreen, which is particularly challenging for MODIS snow detection due to a combination of two factors: i) it is not clear how well canopyintercepted snow is detected by optical sensors [40,41], and ii) the density of the canopy hampers remote snow detection underneath forest [23,42].
Forest-free areas are characterized by a larger mean absolute error during snowmelt season, especially during late melting (i.e., July to September). This occurs because most of the forest is at intermediate elevations (see Table 1), where snowmelt starts around mid-April and totally disappears before summer (see Figure 5c). This means that fSCA estimation in summer at low-intermediate altitudes (i.e., forest sites) is relatively easy since little snow is still present. However, if this portion of the catchment is removed from the estimation of fSCA (like for fSCA of Figure 5b), high-elevated sites play a major role in MAE estimation, since snowpack is still present there, leading to a larger MAE. Moreover, MODIS snow detection during summer months is known to be affected by several uncertainties such as patchy snow cover and clouds/snow misclassification. Indeed, patchy snow cover causes an overestimation of the actual SCA by MODIS, which tends to detect a continuous SCA according to [14,17]. Besides this, summer storms triggered by orographic effects lead to frequent cloudy conditions close to the mountain tops, and this is another potential reason for MODIS fSCA overestimation since clouds might be misclassified as snow [21].
Moreover, by analyzing the snow depth recorded at the 28 snow stations across the catchment, it was found that from July to September, only stations above 2750 m detected snow. A fraction corresponding to the 6% of the total catchment area has equal or higher elevation, which seems to further confirm that the average MODIS-based SCA found over July to September (on average equal to 12.5%) is overestimated.
Concerning the TOPMELT fSCA uncertainty band, Figure 5 shows that it is larger during melt seasons than the accumulation periods. This occurs since the most sensitive parameter is the CMF, controlling most of the SCA dynamics driving the snowmelt rate when temperature exceeds a fixed threshold (see Equation (1)). During the accumulation period, parametric uncertainty is reduced, especially from January to March 2014, when intense snowfall occurred, leading to a nearly complete snow coverage of the catchment. In such cases, it is speculated that parametric uncertainty might have an impact more on the SWE spatial distribution rather than on an aggregated variable such as the fSCA.
In order to investigate whether any spatial pattern affects the agreement between MODIS maps and TOPMELT simulation, observed and simulated SCA are compared using the Overall Accuracy index (OA), following [17]. OA index is estimated pixel-by-pixel over all periods covered by MODIS maps, i.e., from 01 Jan 2011 to 30 Sep 2016, as follows: where: i is the pixel over which OA is computed, Δt is the time window, Nagreements are the times when a pixel i is equally classified as snow-covered or snow-free from both MODIS and TOPMELT and Ntimestep is the number of comparisons between observations and simulations within Δt. When OA is equal to 1, it means perfect agreement throughout the time window, while 0 indicates total disagreement. Moreover, OA analysis is performed aggregating MODIS and TOPMELT snow cover maps by month, in order to investigate the variability between MODIS and simulated SCA depending on the season. OA is computed for each behavioral parameter set, and then the simulation corresponding to the OA median is taken as a reference value for the TOPMELT estimation. Figure 6 shows the monthly pattern of OA in terms of median and interquartile range, depending on different land uses. The interquartile range is meant to quantify the spatial variability of OA for each specific land cover class. July, August and September were not considered because of the little snow coverage across the catchment. Therefore, these months are not reported in Figure 6. An evident seasonal variation of OA can be observed for each land cover class, but for agriculture, which shows very good agreement throughout the year. This happens since agriculture land cover takes place at low elevations, where snow dynamics are not very significant (i.e., snow cover is zero most of the time).
Forested areas during January and February show the lowest OA index among all months and all land cover classes, i.e., the poorest agreement between TOPMELT and MODIS-based SCA ( Figure  6). However, even in such months, a fraction of the forested area shows high OA, resulting in a median value slightly above 0.6 in January and slightly below 0.7 in February. This suggests that the observed bias might not be totally related to the shadowing effect of the forest canopy hampering an optimal snow detection by MODIS. This finding partially agrees with [14,15], which reported satisfactory performances of MODIS even under forest canopy, even though most of the literature clearly highlights a lack of reliability of MODIS-based SCA over forest areas.
During snowmelt season (March to June) the OA estimated over forest areas depicts an increasing trend, showing a distribution of OA in April equivalent to the one resulting on grassland. Then, from November, when the accumulation season begins, OA over forest sites starts again to decrease.
Concerning grassland and bare soil land covers, the OA pattern across the year is similar, with a minimum median value, respectively, in May and in June, when most of the snowmelt occurs. For the bare soil, it happens one month later due to the higher mean elevation (see Table 1), which delays snowmelt peak. However, the minimum OA index over these two land classes is greater than the minimum that resulted over forests, confirming the challenge of providing reliable satellite-derived SCA estimates there.
Regardless of the land cover, the OA index shows large spatial variability as proven by the interquartile range provided in Figure 6. This is mainly related to the elevation range that encompasses each land use class. For instance, the highest values of OA found over forested areas during January and February are likely due to the complete snow coverage that occurred on the forest sites at the highest elevation. In this case, despite the uncertainty related to both TOPMELT and MODIS, reproducing the correct SCA is relatively easy. On the other hand, on forest sites placed at low-intermediate altitudes a more heterogeneous and dynamic snowpack occurs, explaining the drop in OA index. The same mechanism takes place for the other land classes. Figure 6 points out that forest sites are the most critical concerning the accuracy of SCA dynamics, especially during winter months. Therefore, a specific analysis over forest sites from January to March was developed and shown in Figure 7. We select such months among the overall winter period because over forest, they doubtlessly display the worst OA. A drop in OA in forest sites occurs also on December, but in such cases, OA values are similar for forest, grassland or bare soil land cover classes; hence, December is not reported in the following in-depth analysis carried out on forest sites.

A criterion to Identify Critical Areas for Satellite-Derived Snow Detection
Focusing from January to March, the Overall Accuracy was related to the solar illumination conditions using the clear-sky solar radiation estimated with TOPMELT to evaluate the level of solar illumination. Since TOPMELT aggregates the overall spatial distribution of solar radiation in 10 classes (see Section 2.3), Figure 7 shows 10 values of OA accordingly. This procedure aims to identify a possible threshold of solar radiation above which the overall accuracy is still acceptable even under forest canopy. Following [17], OA greater than 0.7 is considered as satisfactory agreement between TOPMELT and MODIS-based SCA. To fit this purpose, the simulated period is divided into two subsets. The first subset (01 Oct 2011 to 30 Sept 2014) was used to quantify such radiation threshold, and the second one (01 Oct 2014 to 30 Sept 2016) allows a validation of the defined threshold. The calibration period allows the identification of a minimum clear-sky solar radiation equal to 17 MJ m −2 d −1 in order to have OA larger than 0.7 on the forest areas (i.e., the first quartile of OA distribution greater than 0.7). The same threshold holds during the validation period, meaning that the minimum solar illumination condition required to have a satisfactory value of OA is not dependent on the period over which it is estimated. Furthermore, Figure 7g-i shows the OA as a function of illumination conditions during the validation period estimated over forest-free sites. OA values corresponding to solar radiation conditions exceeding the identified threshold show negligible difference with OA values found for the same solar radiation levels over forest sites during both calibration and validation sub-periods. Conversely, OA corresponding to low illumination levels (i.e., clear-sky solar radiation lower than 17 MJ m −2 d −1 ) displays generally higher OA than forest areas, especially considering January ( Figure  7a,d,g) and February (Figure 7b,e,h).
In general, the lower the illumination level, the lower the OA, since larger-view zenith angle combined with the presence of trees reduces the viewable gap fraction hindering reliable MODIS acquisition [15,22,24].  Accordingly, Table 4 supports this finding reporting the MAE estimated throughout the available period (01 October 2011 to 30 September 2016) for the different seasons. During the snowaccumulation season, the MAE found applying the proposed criterion is consistently lower than the MAE estimated considering the overall catchment, thus confirming the benefit of removing lowilluminated forested sites from the fSCA comparison. Moreover, it is worth noting that the fSCA MAE found on the accumulation season is 7.5% when the proposed criterion is applied (see Table 4), while the MAE estimated on open sites during the same season is 7.1%. Therefore, applying the proposed criterion leads to an fSCA MAE similar to the one obtained removing the forest sites, but with the advantage of not masking out a larger fraction of the catchment (see Table 5).
Indeed, the main benefit of such a criterion consists in saving part of the overall forested area for the comparison, which represents a potentially significant fraction of the catchment area. For instance, in our case the forest covers 32.17% of the area (see Table 1). The application of the described criterion leads to the identification of the fraction of the forest sites deemed to be acceptable through the year (Table 5). In particular, the period from October to April is shown, since this is the typical snow accumulation and melting period on forest sites, according to the climatology of the region. Likewise, for the period of interest (2011−2016) the snow depth records available across the catchment included in the elevation range where forests are present (i.e. averaged elevation about 1500 m, see Table 1) confirms that snowpack dynamics occurs over such time interval. Despite snow stations are placed on open sites and therefore forest canopy might lead to different snowpack conditions [43], we believe that the seasonality of the snow cycle is the same.  Table 5 shows that in January, a relatively small, yet not negligible, the fraction of the forested area matches the proposed criteria in terms of minimum solar illumination level. Then, this fraction increases up to 71% in April, meaning that MODIS-based SCA estimates can be considered reliable over most of the forest sites during a crucial month for snowmelt in the considered catchment. This effect is expected to be even more important for smaller catchments at high-intermediate altitudes due to the higher fraction of forest land cover.
MAE shown in Table 4 can also be compared with the ones reported in Table 3 in order to understand how TOPMELT fSCA reproduces the satellite-derived fSCA out of the TOPMELT calibration period (i.e., from 01 October 2012 to 30 September 2014; see Table 3). Considering the overall catchment scenario, no remarkable differences are reported, even though during early snowmelt season, the overall MAE (Table 4) was slightly larger than the MAE obtained during the calibration period (9.8% vs. 7.9%). This is expected, since out of the calibration period, the statistical metrics involved in model performance evaluation usually get worse. Nevertheless, the early snowmelt season MAE resulting over the 5-year period is 9.8%, which is still an acceptable value according to the findings of [44,45]. Moreover, the MAE during late snowmelt season is lower during the 5-year period than the one during the calibration period, while during winter, the MAE difference is negligible. These contrasting differences in MAE between calibration (Table 3) and the overall 5year period occurred likely because of either parameter and/or data uncertainties rather than model structure issues.
A further look at MAE of early and late snowmelt seasons reported in Table 4 shows that both MAEs do not show significant differences depending on whether the criterion is applied or not. This happens because, during such months, the area removed following the criterion is limited, since the solar illumination is larger than the threshold of 17 MJ m 2 d −1 over most of the catchment.
Indeed, the proposed criterion contributes to better understanding the reliability of satellitederived snow products during winter months. This might help the assessment of key hydrological variables of snow-dominated catchment, e.g. the maximum snow accumulation. In particular, estimating such variable is often uncertain due to the lack of snow-related ground measures across the catchment [46] and the high spatial variability of snow height [2]. This uncertainty, in turn, affects the correct estimation of snowmelt runoff [47] as well as summer low flows [48]. Thus, conditioning hydrological models using consistent satellite-derived products might be an appealing way to improve hydrological modeling reliability, as proven by [12,49,50] among others.
Finally, it should be noted that the criterion here proposed to identify challenging areas for satellite-derived snow cover acquisition requires neither a large amount of data nor sophisticated implementation, since it based solely on land cover classification and modeled clear-sky solar radiation. Therefore, no particular issue should arise when this procedure is applied in other regions worldwide. However, further investigations are needed when dealing with forest types different from evergreen. For instance, a forest-specific solar radiation threshold might be necessary because of solar radiation attenuation varies depending on different canopy structure.
Furthermore, another potential source of uncertainty in this work is the spatial representativeness of meteorological data, which might be improved using a temperature lapse rate variable in space. This might help in capturing some local phenomena related to the complex topography.
Conversely, model parameter uncertainty was explicitly taken into account through a multi-site model conditioning using GLUE. The multi-site approach is recommended in order to capture the snow dynamics' spatial variability at catchment scale. However, TOPMELT is suitable even for small scale applications, since its spatial resolution can be scaled according to the user's needs. This opens interesting opportunities to take advantage of different sources of reference data, ranging from highresolution satellite data (e.g. LANDSAT or Sentinel-2) to snow cover retrievals of digitaly imagery, which can provide extremely detailed information at the local scale, [51,52].

Conclusions
This paper presents an evaluation of MODIS-based SCA reliability on forest sites, carried out by means of a comparison with TOPMELT, a recently proposed distributed snowpack model. TOPMELT parameters were defined by applying the GLUE procedure through a multi-site comparison of simulated and observed snow depth at snow stations across the catchment. Then, the Overall Accuracy (OA) index was used to investigate land use control on the agreement between simulated and satellite-derived SCA. This analysis pointed out the significant impact of forest sites (mainly evergreen in the studied case), especially during winter months. However, this poor agreement is not uniform within forest areas, but it is strictly related to the solar illumination conditions. A simple criterion was developed in order to identify areas where the combination of low solar illumination level and forest does not guarantee a reliable remote-sensed detection of the snow cover underneath the canopy. This criterion depends on the land use classification and the clear-sky solar radiation distribution solely, both of which were already adopted by several existing snowpack models. In particular, for the examined catchment, a minimum clear-sky solar radiation threshold equal to 17 MJ m −2 d −1 was found, which allows identifying a fraction of the total forested area ranging from 10% in January up to 70% in April, over which MODIS-based SCA is deemed to be reliable despite the presence of the forest canopy.
Finally, the proposed criterion might help with a more consistent use of MODIS retrievals for multi-objective calibration and data assimilation purposes, although it is advisable to perform further testing in different regions, characterized by different climate and forest types. Alongside that, it would be interesting to test the proposed criterion against other satellite-derived products with different spatio-temporal resolution than MODIS, such as LANDSAT or Sentinel-2.