## 4 Discussion  
In this study, three approaches were compared, differing in the type of input data used, based on the same two-stage workflow: the first aimed at filling the gaps of the MODIS snow time series and the second at performing data fusion between MODIS and S2. In both steps, random forests were grown using the satellite MOD10A1 and S2 snow map datasets as training, respectively. Therefore, it must be noted that the training dataset is not entirely free of errors and uncertainties, e.g., due to forests and shaded areas, which have not been explored in this study. The proposed strategy was applied in one region, extending to about 2066 km\({}^{2}\), of the European Alps, and included 7 years of data. As the entire workflow is fully automated and requires only readily available satellite data and a DEM, it allows for applications to other study areas and for several years. This comparison shows that all three approaches are valid, as supported by the high accuracy metrics values, and allow the reconstruction of a daily time series of snow cover at high resolutions, using various types of data: FSC, SCE, and NDSI maps. Overall, we can say that, among them, the A1 approach, which takes as input the FSC MOD10A1  
Figure 10: Variable importance of the regression random forests respectively used in step 1 (a), used to retrieve the gap-filled MODIS now data time series, and in step 2 (b), used to retrieve the daily high resolution snow cover time series. dataset in the first step, and the SCE product in the second, offers the best opportunity, both in terms of accuracy and computational resources. Compared to the other two approaches, A1 is faster, since it does not require a further step to extract the binary snow map from the NDSI map with an optimized threshold and the masking of water bodies, and, moreover, the classification RF is inherently faster and lighter compared to the regression RF performed in step 2 by A2 and A3. When looking at the individual daily maps retrieved, outliers are present, but the entire time series is robust, as confirmed by the results. A future improvement of the algorithm, however, could introduce a probability score or post-processing smoothing based on contiguous days to handle these outliers. To keep the algorithm easily reproducible in other areas, readily available variables were used: orographic features retrieved from a high-resolution DEM, S2 and MODIS data, which are globally and freely available. On the other hand, variables used in other studies, such as Land Use/Land Cover (LULC), distance to water bodies, and variables related to snow transport, such as wind fields and distance to topographic barriers, as in [PERSON] et al. [58], were discarded. This choice was guided by the fact that this work aims at reconstructing long time series of snow; thus, this would require high-resolution time series maps of LULC to be consistent, which are difficult to obtain. Moreover, the LULC also suffers from unpredictable errors, especially if moving backwards in time, so it was deemed appropriate not to use it. Other studies used the LULC as a variable for retrieving snow cover maps [96; 97], but, in those cases, the data resolution and availability were consistent with the analysis scale. At the European level, CORINE Land Cover cartography at 30 m of resolution is available, which, however, is updated every 6 years. In general, the availability of LULC time series at high resolution represents a data gap at a global level. Despite this, the overall accuracy achieved by [PERSON] et al. [58] is slightly higher compared to our results (from + 1 up to + 5% based on the different seasons), while the F1-score is practically the same. Wind speed is an important parameter for the simulation of fine-scale snow distribution [98; 99; 100]; nevertheless, in [PERSON] et al. [58], this variable has a low importance compared to the others. This is probably due to the resolution of the data used by the authors not being consistent with the scale of analysis.