## 2 Materials and Methods
### Workflow: Two-Stage Random Forests  
The implemented workflow, based on random forest (RF) [70], aims to obtain a complete high-resolution time series of snow cover maps by performing data fusion between MODIS and S2 images. The algorithm was developed in R programming language and is based on two steps: the first is aimed at gap-filling the MODIS snow dataset, which is then downscaled and merged with the S2 data in the second step (Figure 2). In the first step, the MOD10A1 snow dataset is used as training variable, while as predictive variables are used topographical features, i.e., elevation, slope and aspect, temporal features, i.e., Day of Year (DOY) and year, and latitude and longitude. The selected topographical variables exert a strong control on snow distribution [58; 71], and the temporal variables are used to exploit the strong interannual consistency of the snow and, thus, reinforce the prediction ability. The use of latitude and longitude also serves as a proxy for local micro-climatic variability. The random forest is constructed using simple  
Figure 2: Algorithm workflow chart. **(a1)** MODIS dataset: in Approaches 1 and 2 is selected the Fractional Snow Cover (FSC), in Approach 3 is selected the Normalized Difference Snow Index (NDSI). **(a2)** MODIS dataset of step 1 after gap-filling. **(b)** Sentinel-2 dataset: in Approach 1 is selected the Snow Cover Extent (SCE) product derived from Sentinel-2 imagery. In Approaches 2 and 3 is se-lected the NDSI derived from the same Sentinel-2 imagery. **(c)** Random forest 1: in all the three approaches is a regression random forest. **(d)** Random forest 2: in Approach 1 is a classification random forest, whose output is directly the binary snow cover, in Approaches 2 and 3 is a regres-sion random forest, which gives as intermediate output the NDSI map, subsequently converted in the final binary snow cover output. random sampling with a size of 25% of each MODIS image, excluding invalid pixels (e.g., cloud cover, missing data, poor light conditions, etc. ), which is, then, used to fill the gaps in each MODIS scene. In the second phase, the same predictor variables are used, including the gap-filled MOD10A1 dataset, and, as training dataset are ingested the aforementioned snow products derived from 200 S2 scenes with cloud cover of less than 10% (Figure 3). Again, the random forest is used to complete missing parts of scenes, due to cloud cover or partial acquisitions, and to retrieve the total scene in days without acquisitions (Figure 4). In the case of cloud cover greater than 10%, the S2 scenes are discarded and reconstructed in their entirety. To find the highest accuracy, we compared three approaches that differed in the selection of data used for random forests training and, consequently, in the type of random forest itself (Table 1):  
* The first approach (A1) uses the FSC data from MODIS as input for a regression RF in the first step, along with the binary Snow Cover Extent (SCE) product from S2 for a classification RF in the second step. * The second approach (A2) uses the MODIS FSC as input for a regression RF in the first step, and the Normalized Difference Snow Index (NDSI) from S2 for the second step. * The third approach (A3) uses the raw NDSI maps from both MODIS and S2 and a regression RF in both the first and second steps, respectively. Figure 3: S2 scenes with a cloud cover threshold below 10% used to train the random forest (**a**) per year and (**b**) per month. The training dataset for the second step RF was built through a stratified random sampling approach, with a sample size of 2.5% of each available scene with a cloud cover lower than 10%, retaining only simultaneously valid pixels both in S2 and MODIS. The stratified random sampling for the construction of the training dataset is designed to uniformly select within each topographic variable considered.