Spatial Downscaling of MODIS Snow Cover Observations Using Sentinel-2 Snow Products

Understanding those processes in which snow dynamics has a significant influence requires long-term and high spatio-temporal resolution observations. While new optical space-borne sensors overcome many previous snow cover monitoring limitations, their short temporal length limits their application in climatological studies. This work describes and evaluates a probabilistic spatial downscaling of MODIS snow cover observations in mountain areas. The approach takes advantage of the already available high spatial resolution Sentinel-2 snow observations to obtain a snow probability occurrence, which is then used to determine the snow-covered areas inside partially snow-covered MODIS pixels. The methodology is supported by one main hypothesis: the snow distribution is strongly controlled by the topographic characteristics and this control has a high interannual persistence. Two approaches are proposed to increase the 500 m resolution MODIS snow cover observations to the 20 m grid resolution of Sentinel-2. The first of these computes the probability inside partially snow-covered MODIS pixels by determining the snow occurrence frequency for the 20 m Sentinel-2 pixels when clear-sky conditions occurred for both platforms. The second approach determines the snow probability occurrence for each Sentinel-2 pixel by computing the number of days in which snow was observed on each grid cell and then dividing it by the total number of clear-sky days per grid cell. The methodology was evaluated in three mountain areas in the Iberian Peninsula from 2015 to 2021. The 20 m resolution snow cover maps derived from the two probabilistic methods provide better results than those obtained with MODIS images downscaled to 20 m with a nearest-neighbor method in the three test sites, but the first provides superior performance. The evaluation showed that mean kappa values were at least 10% better for the two probabilistic methods, improving the scores in one of these sites by 25%. In addition, as the Sentinel-2 dataset becomes longer in time, the probabilistic approaches will become more robust, especially in areas where frequent cloud cover resulted in lower accuracy estimates.


Introduction
The spatio-temporal evolution of snow cover in mountain areas is relevant for understanding runoff production and improving water management in mountain areas [1][2][3]. It is also important for studying plant and animal phenology [4][5][6], soil erosion [7] and glacier mass balance [8,9]. Moreover, from the global energy budget perspective the snow cover has a determinant role since it can cover large areas and has a very high albedo [10,11]. The importance of snow dynamics has motivated the use of a wide variety of observation In this study, a simple methodology is proposed to spatially downscale MODIS snow cover mapping. Through training with Sentinel-2 observations, two approaches for computing probabilities are proposed and evaluated. The methods rely on the high interannual repeatability of snow distribution patterns [48,49], thanks to the strong topographic control of snow distribution at scales below 100-200 m, which remains similar throughout different snow seasons [50][51][52][53]. The two methods exploit the already available database of 20 m spatial resolution Sentinel-2 snow cover maps from the Theia snow collection [18]. Both methods compute the probability of snow occurrence for each Sentinel-2 pixel over particular time periods (seasonal computation), but one of them calculates probabilities in Sentinel-2 pixels independently of each other and the other computes the snow occurrence inside every single MODIS pixel (500 m grid cell size). These probabilities are then used to determine the snow distribution inside partially snow-covered MODIS pixels. The final outputs of these products are 20 m spatial resolution snow distribution maps with a daily frequency, only limited by cloud presence. The two methods were trained and evaluated in three study sites in the Iberian Peninsula with dissimilar characteristics. This evaluation was performed through a leave-one-out validation procedure.

Study Sites
The three study sites are located in the Iberian Peninsula: the Pyrenees, Sierra Nevada and Picos de Europa ( Figure 1). All three exhibit long-lasting snow cover over very complex topography. Nonetheless, their climate and snow dynamics are highly contrasted. The study site located in the Pyrenees encompasses the National Park of Ordesa y Monte Perdido, with elevations between 1060 m a.s.l. (above sea level) and 3350 m a.s.l. This site is characterized by steep U-shaped valleys due to past glacier activity [54,55]. The study area lies in the central Pyrenees between France and Spain ( Figure 1). The Pyrenees are characterized by a high interannual variability of snow accumulation [56], and, although high-elevation areas have a remarkable presence of snow throughout the year, lower zones may undergo repeated accumulation and melting events during the snow season [57].
The Sierra Nevada mountain range is located in south-eastern Spain ( Figure 1). This range hosts the highest summit of the Iberian Peninsula (Mulhacen Peak), which is less than 40 km from the Mediterranean Sea. In spite of the long-lasting presence of snowpack above 2500-3000 m a.s.l., patchy snow distribution in particular periods is frequent even at high elevation [58]. The elevation gradient of this site ranges from 1230 m a.s.l. to 3479 m a.s.l. This study site comprises nearly the entire National Park of Sierra Nevada.
The National Park of Picos de Europa is a mountain massif in the Cantabrian mountain range. It has a moderate elevation, ranging from 530 m a.s.l. to 2650 m a.s.l. Nonetheless the humid climatic conditions of Picos de Europa produce a long-lasting snow cover in the higher areas [59]. Landscape here is very abrupt, being characterized by steep slopes and numerous large cliffs at the highest elevation.

Probabilistic Spatial Resolution Increase
The methods described here exploit two satellite imageries: the Normalized Difference Snow Index (NDSI) of MODIS 6 and C6.1 products [23,60] and the binary snow occurrence maps of the Theia snow collection from Sentinel-2 [18]. Cloudy areas are masked in the NDSI imagery downloaded from the EathData NASA database [60]. Unlike the NDSI of MODIS, which is a measure of the relative magnitude of the characteristic reflectance difference between the visible and shortwave infrared reflectance [22], Theia is a binary database rendering 1 for snow-covered pixels, 0 for snow-free pixels and NoData in case of cloud (or cloud shadow occurrence). From the available Theia imagery, the snow occurrence probabilities were derived as described below. To provide an easier interpretation of the snow-covered area fraction of MODIS pixels, the Fractional Snow Cover (FSC) was derived from the NDSI through a universal equation (Equation (2) in [61]). Two approaches are proposed to determine the snow occurrence probability (Prob 20m ). Later Prob 20m is used to classify partially snow-covered MODIS pixels downscaled to 20 m as snow covered or snow free. One approach computes the snow probability inside each MODIS pixel (hereafter ProbMODIS) and the other is computed for the entire Sentinel-2 tile (hereafter ProbSentinel).
ProbMODIS: This procedure yields a snow probability occurrence (Prob 20m ), based on the snow distribution observed when MODIS pixels were partially snow-covered during a specific training period. The method computes Prob 20m inside each native MODIS pixel (500 m), independently of the rest of the MODIS pixels. For MODIS and Sentinel-2 contemporary dates, the snow occurrence probability was computed for all cloud-free MODIS pixels having completely cloud-free Sentinel-2 pixels inside, and in which the MODIS FSC was between the upper and lower thresholds. Prob 20m was simply the result of summing up the number of times that each Sentinel-2 pixel was snow-covered and dividing by the total number of times that the entire MODIS pixel was partially snow-covered. This probability may also be interpreted as the probability of a Sentinel-2 pixel being snowcovered divided by the probability of a MODIS pixel being partially snow-covered (see Appendix A for further assessment on probabilities).
ProbSentinel: In this case, only the Theia snow collection images for a particular time period were used to determine Prob 20m . This method computes the probability of a particular pixel being snow-covered by simply summing up all days with snow occurrence and dividing this value by the total number of cloud-free observations for this pixel. ProbSentinel probability matches the definition of the Snow Persistence index proposed by [30] for Landsat-8 and Sentinel-2 observations.
The code with a complete example to compute ProbMODIS and to determine the snow occurrence distribution from daily MODIS FSC is provided in the repository https://github.com/ealonsogzl/ProbModis (accessed on 10 October 2021). All cloudcovered MODIS and Sentinel-2 pixels are considered as NoData and thus are not considered for computing probabilities, nor for deriving the spatial downscaled MODIS snow obser-vations. Nonetheless ProbMODIS routines allow the filling of missing Terra values due to cloud presence with Aqua observations of these are available for the same acquisition day, the option used here after in this research.
The final output of the methodology consists of 20 m resolution binary (BinSnow 20m ) snow cover maps (Snow (1)/NoSnow(0) or NoData). First, 500 m MODIS pixels are resampled to the 20 m Sentinel-2 grid with these requirements; if the FSC of an original MODIS pixel exceeds or does not reach a given upper or lower threshold, all BinSnow 20m pixels inside will be respectively considered as snow-covered or snow-free. For FSC between the upper and lower thresholds (determined with a sensitivity test for each test site), the probability of snow occurrence for each 20 m pixel (Prob 20m ) is used to assign pixels as snow-free or snow-covered. Hereby, Prob 20m probabilities inside each 500 m pixel are sorted, and the higher probabilities are then assigned as snow-covered until the number of snow pixels renders the FSC. If two or more Prob 20m have the same probability value, a second sorting criterion is applied. This criterion sorts pixels with the same probability by elevation, from higher to lower, and those with higher elevations are assigned as snowcovered. Hence, in partially snow-covered MODIS pixels, the sum of all 20 m snow-covered pixels (BinSnow 20m =1) divided by the total number of Sentinel-2 pixels inside the original 500 m MODIS raster is identical to the original FSC value. Figure 2 shows the flow chart for generating the spatially downscaled MODIS snow observations.

Evaluation Metrics
The performance of both methods was evaluated with three error metrics: the Kappa coefficient, the percentage of agreement (Agree), and the total snow-covered area error. All evaluations described below have excluded forested areas, which were masked, and were computed for days having concurrent MODIS and Sentinel-2 acquisitions.
The Kappa coefficient [62] is a broadly used coefficient in the remote sensing community [24,[63][64][65][66] and is widely used in a natural resources context [67]. This coefficient evaluates the agreement beyond chance between a reference and the target to be evaluated. The Kappa coefficient is estimated from a confusion error matrix [68]. In a binary classification, such as that evaluated here, the confusion matrix summarizes all possible scenarios with four elements: true positive, false positive, true negative, and false negative, which respectively correspond to a pixel classified as snow when there is snow, classified as snow when there is no snow, classified as no snow when there is no snow, and classified as no snow when there is snow. From the proportion of cases correctly classified (P o , which may be interpreted as the overall accuracy) and those correctly classified by change (P e ; [62]), the Kappa coefficient is derived using the expression; The percentage of agreement (Agree) is computed by dividing the number of times that two variables agree by the total number of samples considered. The computation of Agree is rather simple and although it does not evaluate the agreement by chance, this estimate is a good indication of the similarity between the product and the reference.
Kappa and Agree were obtained with the irr package R (https://cran.r-project.org/ web/packages/irr/irr.pdf accessed on 10 October 2021), omitting missing data and thus avoiding the evaluation of cloud-covered pixels in any products compared. These two indexes allow a good evaluation of the temporal and spatial evolution of the snow-cover mapping performance of ProbMODIS and ProbSentinel.
In order to evaluate the improvement of the spatial resolution increase based on the snow occurrence probability and a MODIS binarization (Snow/NoSnow), these error estimates were also obtained for increased spatial resolution MODIS images (hereafter BinMODIS 20m ). This directly rescales to 20 m the 500 m MODIS values using a nearestneighbor algorithm and then classifies the 20 m pixels as snow-free or snow-covered based on an FSC threshold (note this threshold is different to those used in the probabilistic binarization). Instead of using a threshold from the existing literature, a sensitivity test was conducted to determine which is the FSC value that renders a better snow cover mapping when compared to the Sentinel-2 snow cover map.
As the only 20 m spatial resolution observation of snow distribution is the Sentinel-2 imagery, and this database is also exploited to generate the probabilities with the two methods proposed here, a leave-one-out cross-validation [69] was applied. This method involves the division of the original dataset into N subsets. The model is trained N times with the N − 1 subsets and tested in the remaining N subset. In this case study, the Sentinel-2 imagery is divided into snow seasons (from early September to late July) and thus the evaluation of the ProbMODIS and ProbSentinel is performed with Sentinel-2 images not included in the probability training.
In addition to the Kappa coefficient and Agree, the absolute snow-covered area error was computed to determine if the probabilistic spatial resolution increase of MODIS snow products (either ProbMODIS or ProbSentinel) shows a marked improvement when mapping at the regional scale. The absolute difference between the Sentinel-2 snow-covered area and that obtained with the two methods and BinMODIS 20m (the directly increased resolution MODIS binary classification) was determined.

Results
The spatial distribution of the snow occurrence probabilities (Prob 20m ) accounts for the weight that topography has on the snow distribution in mountain areas. This can be observed for the spatial distribution of ProbMODIS probabilities computed for the entire database of Sentinel-2 imagery in the three study areas (Figure 3). In this way, valleys, ridges and mountain summits can be identified in the snow occurrence probability maps.

Sensitivity Test of FSC Thresholds
The sensitivity test to determine the optimal FSC threshold ( Figure 4) to classify the BinMODIS 20m pixels as snow-free or snow-covered, yielded 0.45 for the three test sites. This sensitivity test comprised all concurrent MODIS and Sentinel-2 images since the later satellite mission became operational. For days with cloud coverage (intersection of MODIS and Sentinel-2 cloud coverage) below 90%, the Kappa value was computed in the cloud-free area. Although a 0.45 FSC threshold is used hereafter to binarize snow presence/absence in the 20 m rescaled MODIS snow cover maps, FSC values ranging from 0.35 to 0.55 render nearly the same results.
The determination of FSC thresholds to consider the 500 m MODIS pixels as snowfree, partially snow-covered or totally snow-covered, in order to then apply ProbMODIS or ProbSentinel probabilities, was assessed through a paired sensitivity test (to compute the upper and lower thresholds). As in the BinMODIS 20m sensitivity test, the Kappa values for both methodologies were obtained in the intersection of cloud-free areas for dates having concurrent Sentinel-2 and MODIS acquisitions with cloud coverage less than 90%. Figure 5 depicts the mean Kappa coefficients for the FSC thresholds in the three test sites and for both methodologies. The sensitivity test shows that for the three test areas different upper and lower thresholds are required. For the Pyrenees, both ProbMODIS and ProbSentinel had lowerupper thresholds of 0.25-0.75 respectively, 0.15-0.85 for Sierra Nevada and 0.15-0.75 for Picos de Europa, respectively. Table 1 summarizes the intervals defined with these thresholds.  Table 1. FSC thresholds to classify the snow cover in the 500 m MODIS pixels. The snow occurrence probabilities are only used in the partially snow-covered pixels. The 20 m pixels inside the NoSnow MODIS pixels area are all assigned as snow-free, and for the AllSnow covered cases, all 20 m pixels are set as snow-covered.  7 show the scatter plot of the two bias estimates, obtained for each day with ProbMODIS/ProbSentinel probabilities versus that obtained with BinMODIS 20m . Nearly all points are above the 1:1 line, showing that the spatial resolution increase based on the snow occurrence probability improves snow-cover mapping. However, for a small fraction of days the methodology fails and produces worse results. These days are generally found at the beginning of the snow season (dark blue) when the snow distribution may not be strongly controlled by the topography. Conversely, late snow season days show a marked improvement in Kappa coefficient values ( Figure 6). Similarly, the percentage of agreement (Figure 7) has a similar behavior and shows an improvement as the snow season advances (dark blue to red colors transition). Again, with a distinct error indicator (in this case the percentage of agreement, Figure 7), it is revealed that the spatial distribution of the snow covered pixels obtained with ProbMODIS or ProbSentinel probabilities better captures the true snow cover distribution. In the three test sites, the two error estimates (Kappa and Agree) for ProbMODIS and ProbSentinel provide better average results than BinMODIS 20m observations, and also lower standard deviations, showing the overall improvement of the methodology ( Table 2). The average improvement obtained for Kappa is at least 10% for the two probabilistic methods, improving the scores obtained in Picos de Europa by 25%.

Spatial Impact When Mapping the Snow-Covered Area Distribution
The improvement of the methods was also been assessed spatially. In this way, for each test area both the Kappa coefficient and the percentage of agreement were computed inside each 20 m pixel including all dates with valid information (cloud-free days). In the three sites, ProbMODIS and ProbSentinel yielded more reliable snow-cover maps than those obtained with BinMODIS 20m . In the Pyrenees (Figure 8), improved results were obtained with the two methodologies when compared to the Kappa coefficients of the BinMODIS 20m maps. ProbMODIS has a better performance than ProbSentinel in the valleys but also in the higher areas. Nonetheless this improvement is minor, and the results were similar for almost the entire area. Figure 8 also shows a major shortcoming of MODIS snow retrieval which was not amended with any of the methods proposed. Snow distribution monitoring on steep slopes with a northerly aspect fails. Thus, in north-facing areas of deep valleys and in steep cirques with a northerly aspect, very low Kappa coefficients were obtained. Results obtained for Picos de Europa ( Figure S1 in the Supplementary Materials) also show a remarkable improvement when comparing the two probabilistic procedures to increase the spatial resolution of snow-cover mapping. The spatial distribution of the Kappa coefficients is again better for ProbMODIS, but ProbSentinel also improves BinMODIS20m. Again, north-facing areas show the inability of the MODIS sensor to retrieve snow-cover distribution. In Sierra Nevada ( Figure S2), the spatial distribution of the Kappa coefficients renders equivalent conclusions to those obtained in the other test areas.
Another evaluation of the methods proposed here is to compute the daily fraction of original MODIS pixels (500 m resolution) having better, same, or worse results with the snow occurrence probabilities (ProbMODIS or ProbSentinel) when compared to BinMODIS 20m . For all dates, the snow-cover mapping obtains better results (not shown) for a large fraction of the 500 m MODIS pixels, with no occurrence of worse performance for any date. This shows that ProbMODIS and ProbSentinel are better able to capture the snow-cover area distribution than MODIS observations directly downscaled, with superior results obtained with the ProbMODIS method. Figure 9 depicts the frequency distribution of the absolute snow-covered area between Sentinel-2, the two probabilistic methods and BinMODIS 20m . Overall, the ProbMODIS and ProbSentinel distributions are displaced towards lower errors, showing almost identical results. Nonetheless the mean improvement (vertical dashed lines in Figure 9) when mapping the total snow-covered area remains marginal. Figure 9. Frequency distribution of the absolute error in total snow-covered area estimation obtained for ProbMODIS (upper panel) and ProbSentinel (bottom panel) when compared to the snow-covered area retrieved with Sentinel-2. The two graphs also depict the BinMODIS 20m frequency distribution. Vertical dashed lines show the corresponding mean snow-covered area error during the study period. These graphs include the results for the three sites.

Discussion
The spatial downscaling of MODIS snow products described here firstly accepted that a universal equation to transform the NDSI to FSC (Equation (2) [61]) is able to provide a good estimation of the snow-covered fraction inside MODIS pixel. Probably, other equations used to transform NDSI values to FSC [61,70], or specifically trained on each test site [71], would yield a more accurate representation of the snow-covered area inside each 500 m MODIS pixel. Nonetheless these would be empirical relations that need specific sensitivity tests to determine the thresholds better capturing the FSC classification. The universal NDSI-to-FSC transformation used here, together with the sensitivity tests for each site, provides a suitable classification of the MODIS retrieved snow-covered fraction, which in turn allows exploitation of the snow occurrence probability products.
The FSC threshold selection was based on the best mean value obtained in the sensitivity tests. These tests showed minor differences for close values around the selected thresholds (FSCthreshold ± 0.1 shows nearly equivalent performance), suggesting that the results may be impacted by small differences when selecting them [24].
For the vast majority of days during the study period, the snow occurrence distribution was better captured using the two probabilistic approaches as better Kappa values were obtained. As the snow season advanced, the improvement was more pronounced. Conversely, in the early snow season, the methodology failed more frequently. This result is related to the topographic control of snow distribution but also to the hysteresis behavior of snow dynamics [72,73]. This is an important drawback for the application of ProbMODIS and ProbSentinel, as the methodology has a lower confidence for observing the snow distribution for the initial snow accumulation events of the season. The snow distribution at the end of the snow season is highly impacted by both processes governing the snow accumulation and processes controlling melting dynamics [51,74], which are summed up as the snow season advances. Thus, at the end of the snow season the snow distribution shows repeated patterns [49,75,76], which allows improved results to be obtained compared with days early in the snow season. Such an outcome of the methodology shows that the main hypothesis accepted here, the snow distribution is strongly controlled by the topographic characteristics and this control has a high interannual persistence, has a stronger validity as the snow season advances. Hereby, it is shown that both, ProbMODIS and ProbSentinel snow occurrence probabilities take advantage of the topographic control of snow distribution [50,52,53,77] to downscale MODIS snow observations. Spatially, the largest improvement of the probabilistic snow occurrence is obtained in lower-elevation areas, and zones with repeated accumulation and ablation cycles during the year. This highlights one of the major ameliorations of the methodology proposed here: improved understanding of snow-cover evolution processes in those areas having a discontinuous and more dynamic snow presence. Higher areas, with longer snow presence, also show better results with ProbMODIS and ProbSentinel when compared to BinMODIS 20m . Nonetheless this improvement is minor as snow is present for longer periods as observed in high-elevation zones. This is evident in the Pyrenees above 1800-2000 m a.s.l., which is considered the elevation from which seasonal snow cover is observed [78]. In contrast, the moderate elevation of Picos de Europa or the melting events that can be expected even at the higher elevations of Sierra Nevada [58], causing a snow cover behavior transition above 2700 m [79,80], means that in these sites the improvement in high-elevation areas is also remarkable. One drawback of MODIS snow-cover retrieval, which is transferred to the final product of the probabilistic spatial resolution increase, is the deficient snow observation in very steep north-facing areas as these are affected by the amount of radiation received by a pixel [44,81].
As the Sentinel-2 mission has been available since late 2015 and the revisit time is 5 days, the evaluation of the methodology described here for days without clouds (or with minor cloud extent) would have led to very few evaluation dates [82]. To avoid this limitation, the evaluation strategy is focused on obtaining the bias estimates in cloud-free areas for both acquisitions, MODIS and Sentinel-2, with a 90% threshold of the cloudcovered area intersection. For some dates, this may have an impacted on the evaluation, as the cloud-free area may be in areas that are completely snow-covered or without snow presence. Moreover, clouds' occurrence in mountain areas is highly impacted by the interaction between rugged terrain and the dominant wind direction [83], producing cloud patterns connected with the study area's topography [84]. This irregular occurrence of clouds may have impacted the spatial evaluation of the probabilistic snow-mapping products, as some locations may have a much longer track record than others. Limitations from both cloud-cover occurrence and dissimilar lengths of snow occurrence in different zones of the study areas will become less important as Theia Sentinel-2 imagery continues to acquire images. Longer observation datasets will allow more reliable probabilities of snow occurrence to be generated, as more cases will be included in the training. For instance, Sentinel-2 pixels with the same Prob 20m values which now are sorted by elevation would be rarer as more snow seasons are available. Thus, it is expected that as long as the Sentinel-2 mission continues to acquire observations, the snow occurrence probabilities will continue to improve their robustness, making the elevation sorting of pixels with the same probability unnecessary.
Overall, the snow distribution observations obtained with the methodology described here are limited by the inherent limitations of optical satellite sensor limitations [17,85]. Hereby, the observation of snow distribution in forested areas or restrictions derived from cloud cover presence are shortcomings that will be intrinsic to the observations derived from this spatial downscaling.
Of the two methods described and evaluated, ProbMODIS showed better results than ProbSentinel. Nonetheless, this difference is minor when compared with the improvement obtained from BinMODIS 20m . This shows that accounting for the snow occurrence inside the 500 m MODIS pixels (as ProbMODIS does) is more precise than accounting for the snow occurrence in Sentinel-2 pixels independently of each other. Thus, use of the ProbMODIS algorithm is recommended to retrieve the snow probability occurrence and then binarize the FSC should be binarized according to the thresholds of the sensitivity tests performed.
The probabilistic snow occurrence approach described here is only able to improve the snow-mapping resolution in partially snow-covered MODIS pixels. Because not all pixels, and sometimes a minor fraction of them, are partially snow-covered (when compared to the entire snow-covered area), neither method has a major impact in quantifying the total snow-covered area. This lesser extent of partially snow-covered pixels also explains the fact that for the vast majority of days, similar results were obtained for a larger portion of the 500 m MODIS. However, the improvement obtained in transition zones can have remarkable importance for monitoring snow hydrology as snow-dominated zones lead to rain-dominated zones [86].
Other spatial downscaling methods of MODIS snow-cover products have also provided reliable and accurate snow maps [45], with even higher spatial resolution than that obtained in this study [47]. These authors were able to downscale MODIS snow-cover information to very high spatial resolution (3 m) with very good results. Nonetheless this latter method required both a Digital Elevation Model (DEM) and a training/validation dataset with the same spatial resolution as that of the final product, and a specific sensitivity test to determine the weights of the terrain parameters and their search distances. The main benefits of the method proposed here when compared to previous works [44], [45,47], is that it only relies on space borne satellite observations of the snow cover without any further requirement on the computation of the spectral reflectance or the calculation of further topographic variables. Hereby, as the spatial downscaling of MODIS snow cover observations with ProbMODIS and ProbSentinel, is rather simple and moreover it exploits already available snow cover observations (Theia snow collection and the NDSI from the NASA Earth Observation database), it has an easy and straight forward use. As the methodology obtains reliable and accurate snow cover patterns its simplicity and low computational requirements are the main benefits of its use.
Although the methodology presented requires probabilistic training with the already available database of Sentinel-2 observations, the approach is suitable for a wide variety of applications requiring high spatial resolution information on snow occurrence before the Sentinel-2 mission started. For instance, studies aiming to understand the ecological cycles of mountain pastures require detailed observations of snow disappearance over a longer time period [4,87]. Similarly, soil erosion in steep mountain areas is tightly related to snow gliding phenomena [7], which may occur at slope scales but also at landscape scales [88] and thus can also benefit from long-lasting and high spatial resolution snow-cover observations. The authors of [30] described and exploited a Snow Persistence index to determine longlasting snow-covered areas (due to wind redistribution and avalanche action) and evaluate snow redistribution models. This latter index was derived from Sentinel-2 and Landsat-8 imagery, in the same way that ProbSentinel probability was computed here, showing that the spatial downscaling of MODIS snow cover observations is valuable for evaluation of high spatial resolution snow cover models. Another research topic that would also take advantage of high spatial resolution snow maps is the monitoring of very small glaciers [89], as these observations would allow an understanding of glacier dynamics in areas under threat from climate change [90]. Another subject that can benefit from this MODIS spatial downscaling is the study of the habitat of some animal species, as snow affects the distribution of forage species [91]. Lastly, the seasonal snow cover at low elevations in temperate mountain regions is predicted to cause important snowpack changes [92]. This demonstrates the importance of high spatial resolution snow monitoring with a long-lasting database [40] as this would allow changes occurring in mountain areas to be understood.

Conclusions
The methodology detailed here demonstrated that it is possible to spatially downscale MODIS snow cover observations by simply determining the snow occurrence probability from a database of higher spatial resolution, in this case study Theia snow distribution maps from Sentinel-2. The amelioration of the snow cover mapping is based on the strong topographic control of the snow distribution that characterizes mountain areas. Using several snow seasons to generate a snow occurrence probability for the original 500 m MODIS pixels (ProbMODIS) or an entire Sentinel-2 tile (ProbSentinel), the methodology is able to improve the snow cover mapping in partially snow-covered MODIS pixels generating snow occurrence maps at a 20 m spatial resolution. Nonetheless the snow occurrence probability computed for each MODIS pixel obtained more reliable snow cover maps and is thus recommended. With the already available Sentinel-2 imagery it is possible to train snow probability occurrences for particular test areas, which after finding the FSC thresholds to classify snow pixels (snow-free, snow-covered or partially snowcovered) allows the generation of snow distribution maps at 20 m resolution from MODIS observations. Nonetheless, these snow observations may be impacted by the inherent limitations of optical satellite retrievals, and will have a lower confidence on the observed snow distribution in the early snow season. The simple computation of this method allows an easy application for a wide variety of research topics requiring a long-lasting database of high spatiotemporal resolution observations of snow distribution in mountain areas. Future work aims to derive ProbMODIS probabilities over the Iberian Peninsula mountain areas and then obtain a high spatial resolution database of snow distribution which may be of interest in other disciplines. Funding: This research was funded by HIDROIBERNIEVE: CGL2017-82216-R, the Grant Juan de la Cierva Incorporación IJC2018-036260-I of the Spanish Ministry of Science and Innovation, the ANR TOP project, grant ANR-20-CE32-0002 of the French Agence Nationale de la Recherche.
Acknowledgments: This work was supported by the Spanish Ministry of Economy and Competitiveness project "HIDROIBERNIEVE: CGL2017-82216-R", and the "Organismo autónomo de Parques Nacionales" with the project "Cartografía de alta resolución especial del manto de nieve y su variabilidad reciente en los PPNN de montaña y los impactos del cambio climático para el horizonte 2050". J. Revuelto is supported by the Grant IJC2018-036260-I. S. Gascoin is supported by the ANR TOP project. The authors thanks the three anonymous reviewers and Norah Wang who edited this article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Probability Assessment on ProbMODIS Product
ProbMODIS probabilities may also be interpreted as the probability of a Sentinel-2 pixel being snow-covered divided by the probability of a MODIS pixel being partially snow-covered. This interpretation can be justified with the Bayes theorem as follow. Let us assume that P(S2) is the probability of an S2 pixel being snow-covered (number of times in which snow is observed divided by the number of times there is an observation without clouds) and P(M) is the probability of a 20 m MODIS pixel being partially snowcovered (nearest-neighbor downscaling). ProbMODIS is defined as the number of times that a Sentinel-2 pixel is snow-covered divided by the total number of times that the entire MODIS pixel is partially snow-covered. This ProbMODIS can be expressed as Although for a nearly negligible fraction of snow-covered Sentinel-2 pixels, there are some occurrences of MODIS pixels with a FSC of 0, it can be assumed that The MODIS-Sentinel-2 observation discrepancies show that in some rare cases there can be acquisition errors, as there is snow observed with Sentinel-2 that is not observed with MODIS for the same day under clear-sky conditions. Assuming the truth of Equation (3), ProbMODIS can also be interpreted as: Equation (4) shows that justifies the probability of a Sentinel-2 pixel being snowcovered divided by the probability of a MODIS pixel being partially snow-covered can also define ProbMODIS.