Bias correction and downscaling of snow cover fraction projections from regional climate models using remote sensing for the European Alps

Mountain seasonal snow cover is undergoing major changes due to global climate change. Assessments of future snow cover usually rely on physical based models, and often include post-processed meteorology. Alternatively, here we propose a 10 direct statistical adjustment of snow cover fraction from regional climate models by using long-term remote sensing observations. We compared different bias correction routines (delta change, quantile mapping, and quantile delta mapping) and explore a downscaling based on historical observations for the Greater Alpine Region in Europe. All bias correction methods adjust for systematic biases, for example due to topographic smoothing, and reduce model spread in future projections. Averaged over the study region and whole year, snow cover fraction decreases from 12.5 % in 2000-2020 to 15 10.4 (8.9, 11.5; model spread) % in 2071-2100 under RCP2.6, and 6.4 (4.1, 7.8) % under RCP8.5. In addition, changes strongly depended on season and altitude. The comparison of the statistical downscaling to a high-resolution physical based model yields similar results for the altitude range covered by the climate models, but different altitudinal gradients of change above and below. We found trend-preserving bias correction methods (delta change, quantile delta mapping) more plausible for snow cover fraction than quantile mapping. Downscaling showed potential but requires further research. Since climate 20 model and remote sensing observations are available globally, the proposed methods are potentially widely applicable, but are limited to snow cover fraction only.


Introduction
Mountain regions store large amounts of precipitation in form of snow and ice, which provide essential water supply for downstream regions, affecting an estimated quarter of humanity (Immerzeel et al., 2020). Global warming resulted in 25 significant changes of the cryosphere with melting glaciers and shifts in the timing and abundance of snow (Huss et al., 2017), which already affected the hydrological cycle (Morán-Tejeda et al., 2014) and will continue to do so in the future (Hanzer et al., 2018). These changes imply consequences in water supplies for domestic use, hydropower, and agriculture.
Seasonal snow cover responds rapidly to climate variability and change, in contrast to glaciers, which are out of balance with https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.
The simplest form of bias correction is the delta-change (DC) approach, where the mean climate change signal (e.g., in temperature) is superimposed on the observation series. However, DC cannot reflect any change in the future distribution of the considered variable. The most widely used approach for bias correction is quantile mapping (QM), which can 65 simultaneously perform downscaling, too. QM matches observed and modelled distributions and the non-parametric variant performs better in reproducing observed climatology than parametric versions (Gudmundsson et al., 2012). An extension is quantile delta mapping (QDM), which is a trend-preserving approach to QM (Cannon et al., 2015), since QM has been shown to modify trends in a few cases (Maurer and Pierce, 2014). The flexibility, performance, and ease-to-use has made QM or QDM a standard approach for national climate change assessments, see, for example, Switzerland (CH2018, 2018 or 70 Germany (Krähenmann et al., 2021).
For assessing future changes in snow cover based on climate model scenarios, two methods are mainly employed. The first is to use downscaled and bias corrected meteorological forcing from climate models to drive dedicated snow or hydrological models (DeBeer et al., 2021;Hanzer et al., 2018). The second is to use directly snow cover output from climate models (see above). However, the availability of long-term high-resolution satellite imagery has enabled a third option: To use remote 75 sensing for bias correction and downscaling of RCM snow cover. To the best of our knowledge, this has not yet been performed. We restrict the study to snow cover fraction only, and not snow depth or SWE, because it is the only snow variable available at high resolution and with high accuracy in mountain terrain, in addition to the method thus being potentially applicable globally.
The aims of this study are to bias correct and downscale snow cover fraction from RCMs using remote sensing observations 80 for the European Alps, and to compare this, for a limited area, to the use of a dedicated snow model forced by downscaled RCM meteorology. The motivation is that the snow cover fraction biases of RCMs were shown to be mainly caused by orography and temperature, partly also precipitation, mismatches (Matiu et al., 2020b), which are systematic biases that can be statistically corrected, while future change estimates should be consistent. While bias correction cannot add information beyond what is contained in the RCM, it can reduce model spread and it can make information on future projections more 85 meaningful: absolute values are better than change estimates both for climate change information and for impact assessments, which often depend on absolute thresholds. By exploiting the morphological dependence of snow cover on topography, downscaling can improve the local spatial representation of RCM snow cover fraction. Finally, the comparison between downscaling RCMs and using a snow model shall highlight benefits and limitations of both approaches.
The study combines the proof-of-concept of applying bias correction and downscaling to snow cover fraction with its 90 application to assess future scenarios of snow cover fraction over the European Alps. The remainder of the paper is as follows. Section 2 introduces the study region and data sets. Section 3 explains the methods used for bias correction and downscaling. Section 4 presents results and discussion, and Section 5 the conclusion.

Data
The study region is roughly what is called the Greater Alpine Region in Europe, which contains the European Alps and 95 spans, in this study, approximately 43-48.5° N and 5-17° E (Fig. 1). The large-scale climatic setting includes influences from the Atlantic Ocean, the Mediterranean Sea, and the European continent. The region is characterized by complex topography with strong elevational gradients. The comparison of statistical downscaling to a dedicated snow model is performed for a small subset, the Ötztal Alps region in Austria (1850 km 2 , 862-3770 m a.s.l., Fig. 1b and 1d); see Hanzer et al. (2018) for a detailed description of the region. 100 https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.

Observed snow cover fraction from remote sensing
As for remote sensing observations, we relied on MODIS (Moderate Resolution Imaging Spectroradiometer), because it 110 offers the best tradeoff between temporal availability (two decades, daily) and spatial resolution (250 m nominal) to perform a downscaling -in contrast to coarser products such as based on AVHRR (Advanced Very-High-Resolution Radiometer), which have a longer period into the past (starting in the 1980s), but are of coarse resolution and less quality for complex mountain terrain than higher resolution sensors. A cloud filtered product was used (Matiu et al., 2020a), which is based on the snow maps developed in Notarnicola et al. (2013a). The processing included temporal and elevational (spatial) filters to 115 remove nearly all cloud coverage. Compared to Matiu et al. (2020a), here we used only Terra acquisitions, thus extending the temporal extent back to 2000. Otherwise, the processing is the same, only the first step (merging of Terra and Aqua) was omitted. Consequently, daily cloud-free binary (snow/land) snow cover maps were available for the complete period 2000-02-24 to 2020-08-23 at 250 m nominal resolution (232 m effectively) in Lambert azimuthal equal-area projection for the domain denoted in Fig. 1. 120 The high-resolution binary snow maps were aggregated into low-resolution snow cover fraction maps that match the RCM resolution of 0.11°, which is approximately 8.6 by 12.2 km for the Alps. Each low-resolution pixel then contained 1961 (37*53) high-resolution pixels.

Snow cover fraction from regional climate models
The EURO-CORDEX ensemble consists of 11 RCMs driven by 8 GCMs from CMIP5 (Coupled Model Intercomparison 125 Project Phase 5); see Coppola et al. (2021) for more information on the general ensemble setup. However, not all models provide all variables and/or all emission scenarios. For this study, we used the RCP2.6 (representative concentration pathway) and RCP8.5 scenarios, where RCP2.6 is likely to keep global warming below 2° C until 2100, while RCP8.5 corresponds to approximately 5° C global warming. Regarding snow cover fraction (SNC), the available ensemble for this study included, for RCP8.5, 6 RCMs driven by 6 GCMs with a total of 29 simulations, and, for RCP2.6, 4 RCMs driven by 5 130 GCMs with a total of 8 simulations (Table S1). The list of used RCMs is CLMcom-CCLM4-8-17, CLMcom-ETH-COSMO-crCLIM-v1-1, CNRM-ALADIN63, IPSL-WRF381P, KNMI-RACMO22E, and SMHI-RCA4. Even though DMI-HIRHAM5 also provides SNC, we excluded it because the SNC values over the Alps are unrealistically low, although snow depths match very good (Matiu et al., 2020b).
The RCM SNC grids were reprojected onto the low-resolution MODIS maps using nearest neighbor resampling in order to 135 have a one-to-one correspondence of pixels in the spatial domain. Nearest neighbor was favored over other resampling methods, such as bilinear, because it preserves the two-sided bounded nature of SNC, which goes from 0 to 1, and thus keeps the same limits, while bilinear resampling can introduce lower maxima or higher minima. Some models display snow accumulation issues (Terzago et al., 2017;Matiu et al., 2020b;EURO-CORDEX Errata, 2021) and affected grid cells were removed based on thresholds on snow water equivalent or snow depth, as described in Matiu 140 https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.
(2020b). These were mostly grid cells with the highest elevation and the number of affected cells was between 0 and 233, depending on model, out of a total of approximately 5000 land grid cells in the study domain (see Figure S1 and Table S2 for location and number of affected cells by RCM). Spatial averages and ensemble means are based on the common subset of pixels available to all models. KNMI-RACMO22E was strongly affected by snow accumulation, and thus we removed it from results that show spatial/ensemble means, since otherwise too many pixels would have to be removed. Consequently, 145 for calculating ensemble means, 23 simulations were available for RCP8.5, but only 4 for RCP2.6.

Snow cover fraction from AMUNDSEN
The fully distributed snow and hydroclimatological model AMUNDSEN (Strasser, 2008), now available as openAMUNDSEN in Python (Warscher et al., 2021), has been previously applied to study the future snow and ice evolution 150 in the Ötztal Alps region in Austria (Hanzer et al., 2018); see Fig. 1b and 1d for the area. AMUNDSEN dynamically resolves the mass and energy balance of snow and ice and has been driven by future meteorology based on EURO-CORDEX RCMs, which includes a subset of the same RCMs mentioned above. The input meteorology has been bias corrected and downscaled using QM to point scale, and further temporally disaggregated and spatially distributed to provide 3 hourly forcing at 100 m horizontal resolution for the whole catchment. For more details, see Hanzer et al. (2018). The modelled 155 snow water equivalents were converted into snow cover fraction using a threshold of 5 mm. We also evaluated a threshold of 15 mm but found differences to be negligible.
For the comparison, we decided to focus on ensemble means, and not compare results between individual GCM-RCM pairs directly. Only few RCMs overlap between Hanzer et al. (2018) and this study. In addition, the QM that has been applied in Hanzer et al. (2018) removes single GCM-RCM biases for the past period and thus renders comparison less meaningful. 160 Finally, the ensemble size is similar, since Hanzer et al. (2018) used 14 GCM-RCM for RCP8.5 and 3 GCM-RCM pairs for RCP2.6, compared to 19 and 3 here, respectively (CNRM-ALADIN63 was removed, because half of the Ötztal Alps area was affected by snow accumulation).

Bias correction of snow cover fraction
We compared four different bias correction methods, which are routinely applied for temperature and precipitation series, in their applicability for bias correcting snow cover fraction: delta change (DC), quantile mapping (QM, Gudmundsson et al., 2012), quantile delta mapping (QDM, Cannon et al., 2015), and multivariate QDM (Cannon, 2018). In all cases, the past 170 refers to period of available MODIS observations, which is 2000 to 2020. For the climate model runs, the historical period, https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.
which goes from 1950/70 to 2005, was merged to the RCP scenario run, which covers 2006 to 2100, in order to have the same common period for the past as available from MODIS (2000 to 2020). Thus, for applying the bias correction, each scenario had its own past time series (i.e., for each RCP scenario), while usually the calibration is performed on the same historical run for all scenarios, usually 1970-2000. However, this was not possible here, since the overlap between historical 175 period and MODIS is only 5 years, which is too little to derive robust distributional estimates of the snow cover climatology.
As future period, we considered 2071 to 2100. The bias correction was applied for the same spatial (0.11°) and temporal scale (daily).
For the DC approach, we calculated the multiplicative change ratios between the past and future from RCMs and applied it to the observations from MODIS. This was done separately for each grid cell and each month. For QM and QDM, we 180 employed the standard routines using empirically derived distribution functions (Gudmundsson et al., 2012;Cannon et al., 2015), as available in the R packages qmap and MBCn, again month-by-month and pixel-by-pixel. For QDM, multiplicative change ratios were used. The multivariate QDM was applied in the spatial domain, thus the multivariate component was to account for the spatial correlation in SNC between pixels. However, we found results to be almost identical to standard univariate QDM, and we do not show it further in results. We assume this similarity to be caused by the high spatial 185 correlation in SNC, and the fact that this spatial correlation is similar in both model and observed series.
Because SNC is bounded not only at the minimum 0 but also at the maximum of 1, the standard QM and QDM algorithms had to be modified. We applied the trace condition, which sets all values below a threshold (here: 0.001, also called trace value) to exact zeros, also to the maximum of 1. In addition, the distribution of RCM SNC contained many exact zeros and ones in comparison to observed SNC, which was more regularly distributed across the [0, 1] interval, caused by the sub-grid 190 variability in observations. This caused problems in estimating and matching the modelled and observed quantiles. To alleviate this issue, we added a random component to all SNC values near 0 and 1, where near means half of the trace value.
The random values were randomly sampled from a uniform distribution with minimum the machine epsilon (the lowest value without rounding issue in floating point arithmetic) and maximum half the trace value, i.e., effectively from the [0, 0.0005] interval. This random component helps in matching quantiles (and thus distributions) but breaks the temporal 195 consistency in the bias corrected SNC time series.

Downscaling of snow cover fraction to binary snow 200
The proposed downscaling approach converts low-resolution snow cover fraction (SNC) from RCMs into high-resolution binary snow cover (snow/land), from which we extracted monthly and annual snow cover duration (SCD). The downscaling is based on the morphological dependence of snow cover (Premier et al., 2021). It uses a conditional probability approach (Dong and Menzel, 2016) to define the relationship between snow cover fraction of a low-resolution pixel and the https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.
probability of a high-resolution pixel being snow-covered or snow-free. We used the 20 years of daily MODIS snow maps to 205 calculate the conditional probability that a high-resolution pixel is snow-covered depending on the SNC of the lowresolution MODIS pixel, which has been aggregated from the high-resolution maps. For this, we summed up the occurrence of snow per high-resolution pixel for a sequence of low-resolution SNC bins of width 0.05 with additional bins at the minimum and maximum to catch nearly exact zeros and ones (i.e., 0, 0.001, 0.05, 0.1, … 0.95, 0.999, 1). See Figure S2 for estimated probabilities by SNC bins for an example area. Probabilities were only estimated, if more than 30 observations 210 were available per low-resolution pixel and SNC bin. Low-resolution snow cover fraction threshold at which the probability that the high-resolution pixel is snow covered exceeds 0.5 (SNCp50). (c) Same as (b) but after imputation based on similarity of probability curves (as well as similarity of low-resolution sub-grid topography). (d) Empirically estimated probability curves that high-resolution pixels are snow covered depending on the low-resolution snow cover fraction. Labels (1) to (5) denote five example high-resolution pixels in panels (a-d).

220
The downscaling is then performed by estimating the SNC value for which the probability of a pixel being snow covered is higher than 0.5, called SNCp50 from now on. For this, the SNC threshold at which the probability exceeds 0.5 was estimated from the empirical relationship by linear approximation. An example is shown in Fig. 2, which shows the expected negative relationship between SNCp50 and elevation, which implies that, as SNC increases from 0 to 1, the snow cover is more likely 225 https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.
to be found going from high to low elevations. But while elevation is the main influence, SNCp50 can vary considerably for similar elevations (e.g., points (1), (3), and (4) in Fig. 2 are approximately the same elevation) due to local terrain factors.
In order to have a unique solution for SNCp50, the relationship between low-resolution SNC and probability of highresolution snow needs to be non-strictly monotonic, which is a physically valid assumption for any given high-resolution pixel. However, because of noise and errors in the MODIS time series, this was not always the case, so we selected the 230 longest non-strictly increasing subsequence to calculate SNCp50. To have robust estimates of this threshold, we removed nearly exact zeros and ones, thus requiring some points that identify the shape of the curve (Fig. 2d). For 36% of pixels, the linear approximation failed to estimate SNCp50, because either no empirical estimates were available (except for ones and zeros) or all were above/below 0.5 (see, e.g., point (3) in Fig. 2d). For these pixels, SNCp50 was imputed in two steps.
First, we applied a similar pixel approach (Li et al., 2020) to select another reference pixel with available SNCp50 that is 235 similar with respect to the sub-grid topography of the encompassing low-resolution pixel and to the high-resolution probability curves. For this, we identified similarity between low-resolution pixels, by calculating the Wasserstein distance (also called earth-mover distance) between the high-resolution elevations in all pairs of low-resolution pixels. The Wasserstein distance is a distance metric applicable to distributions, where no one-to-one correspondence can be expected.
For each high-resolution pixel with missing SNCp50, we selected the 50 nearest low-resolution pixels (including the low-240 resolution pixel with missing high-resolution SNCp50). Then, we calculated the mean absolute error (MAE) between SNC probability curves for pixels that deviate at most 150 m of the missing pixel elevation and have at least 5 values to compare SNC probability curves. The SNCp50 from the pixel with minimal MAE was used to fill the gap. After this step, 29% remained missing, most located above 3000 or below 500 m.
The second step involved a simpler elevation filter, and no comparison of probability curves. Again, we selected the 50 most 245 similar low-resolution pixels, and from these up to 100 high-resolution pixels with at most 150 m elevation difference to the gap pixel. The average SNCp50 from these maximum 100 pixels was then used to fill the gap. This step left <0.001% pixels missing.
A rough estimate of the downscaling performance was calculated by applying the downscaling to the QM or QDM bias corrected past model runs, which, by definition, have a matching empirical distribution in SNC to the low-resolution MODIS 250 SNC. Thus, the difference between downscaled average annual snow cover duration (SCD) and observed high-resolution MODIS SCD is an indicator of the downscaling error. The mean downscaling bias was -3.7 d and the MAE 5.7 d. In addition, there was an elevation dependence of the bias ( Figure S3). A negative bias was found for elevations below 1000 m, almost no average bias between 1000 and 3000 m, and positive bias above 3000 m. Glacierized surfaces exhibited strong positive bias, except if SNCp50 was imputed by the second elevational step. We excluded glacierized pixels with more 10% 255 glacierized area from the further analysis of the downscaling, because of this systematic bias in addition to the difficulties of distinguishing snow and ice with MODIS (Fugazza et al., 2021). In addition, they had a strong overlap with the already removed low-resolution RCM pixels with snow accumulation (Sec. 2.1). Glacier extents were extracted from the Randolph Glacier Inventory 6.0 (RGI Consortium, 2017). The downscaling was applied to both QM and QDM bias corrected low-resolution snow cover fraction. However, in the 260 results we put emphasis on QDM, because the trend preserving attribute is preferred for ratio variables, such as precipitation and SNC, as compared to interval variables, such as temperature (Cannon et al., 2015).

Bias correction and future changes in snow cover fraction
The RCMs reproduced overall seasonal and large-scale spatial patterns of past snow cover fraction well (cf. RAW in Fig. 3 and 4). For instance, Winter (December-February) snow cover fraction spatial patterns agreed not only for the high-elevation Alpine region, but also for lower elevation mountains, such as the northern Alpine foreland, Dinarides, or the Northern Apennines. However, because of the coarse resolution and smoothed model orography, the RCMs did not capture the fine 270 scale complex patterns found in the Alps (Fig. 4). Additionally, they revealed some spread in under-and overestimating snow cover fraction, which depended both on RCM and GCM (Fig. 3), as has also been shown previously (Terzago et al., 2017;Matiu et al., 2020b). This model bias depended stronger on the RCM and only secondly on the driving GCM.
Applying DC, QM and QDM bias correction to past RCM output enforced it to match the distribution of observed SNC and consequently also reduced the model spread for the future (Fig. 3). In addition, it introduced the fine-scale spatial patterns 275 into the smoothed model output (Fig. 4). QM and QDM, by definition, resulted in the same patterns for the past. Bias corrected future estimates were similar for the two trend preserving approaches DC and QDM, which themselves differed substantially from QM. For example, QM showed less reduction in SNC under the RCP8.5 scenario for spring (March-May) than DC and QDM ( Fig. 3 and 4).   Projected changes until the end of the century depended strongly on elevation, and the strongest reductions in winter SNC were observed between 600 and 2000 m and in spring above 1000 m ( Figure S4 and S5). Under RCP2.6, winter SNC decreases approximately 7 pp at 1000 m elevation, while above 2000 m it's less than 2 pp and the model uncertainty includes no change (Table S3) (Table S4). In addition, a high model uncertainty in projected changes of spring SNC above 2000 m under 315 RCP8.5 was observed: The model spread ranged from almost no change to an approximate halving of SNC ( Figure S4 and Table S4).
This model spread in spring SNC under RCP8.5 is likely caused by the snow schemes in the climate model's land surface schemes in combination with the projected temperature and precipitation changes, which directly affect SWE, and, since SNC is parametrized on SWE, also SNC. Higher uncertainties are expected in spring because potential errors accumulate 320 over the snow season. However, a detailed discussion on snow model processes and uncertainties is beyond the scope of this study, and better addressed in dedicated projects, such as ESM-SnowMIP (Krinner et al., 2018).
An in-depth view of the bias correction results at a single grid cell highlights the main differences between raw climate model SNC and observations as well as between QM and QDM ( Figure S6). RCMs have more saturated SNC at both 0 (snow free) and 1 (snow covered) and thus display a more of a fully snow free and fully snow-covered grid cell over time as 325 compared to MODIS. This is likely caused by sub-grid variability, which is prominent in MODIS, since it's based on 250 m information. The trend-preserving attribute of QDM keeps the distribution of SNC identical between past and future when raw model SNC does not change, e.g., for the fraction of time, where the grid cell is fully snow covered (Fig. S6). In the same situation, QM shows reductions in SNC. Additionally, QM has spurious breaks caused by applying the method monthby-month, but not QDM. While these breaks could be alleviated by applying the bias correction with a moving window 330 approach (e.g. 3-6 months) or using the whole year, QM still suffers from artificial modification of trends, as has been partly shown before for precipitation (Maurer and Pierce, 2014;Maraun, 2013). Consequently, it should be treated with caution for snow cover fraction, too. For the downscaling below, we thus only used results from QDM.
The European Alps have a prominent north-south climatic divide (Auer et al., 2007), which manifests itself in snow cover 360 duration, too. Taking anomalies of SCD by elevation shows, on average, higher SCD north of the main ridge and lower SCD south (Fig. S9a). These patterns were reproduced in RCMs, too, and changed in the future period: Comparing RCP2.6 to RCP8.5, the north-south gradient in SCD became less strong for lower elevations and more pronounced for higher elevations (Fig. S9). In addition, a stronger relative decline in SCD was observed south and west of the Alps compared to north and east ( Fig. 5c) under RCP2.6. An analysis of station snow depth and SCD trends over the last five decades in the Alps similarly 365 showed stronger declines south than north (Matiu et al., 2021). Consequently, this trend might continue in the future given the findings in this study.
The average downscaling uncertainty (bias and error, see Sec. 3.2, Fig. S3) was much lower than model spread and future change estimates. Thus, the largest part of uncertainty of future projections was less because of the downscaling method, but more caused by the spread in GCM forcing together with RCM snow schemes. Since there is no single "best" climate model 370 (Vautard et al., n.d.) and no single best snow model (Etchevers et al., 2004;Rutter et al., 2009;Menard et al., 2021), we conclude it is safe to take model spread as representative of model uncertainty for future projections. However, the negative bias from downscaling caused future estimates to be also negatively biased as compared to estimates from bias correction only (Fig. S10). For example, the average difference over all low-resolution grid cells in future annual SNC under RCP8.5 between the bias corrected ensemble mean and an upscaled SNC after downscaling was -1.4 pp, which corresponds to -4.9 d. 375 The employed statistical downscaling method extrapolates information beyond the elevation coverage of the RCMs. At 0.11°, and not considering the pixels with snow accumulation, the highest grid cell from low-resolution RCMs was at approximately 3000 m, which contained single high-resolution pixels with elevations up to 4105 m. The downscaled estimates above 3000 m thus should be treated with caution, even though the observed stronger reductions at elevations above 3000 m from downscaling are similar to the results of the simulations from a high-resolution RCM, which explicitly 380 resolved elevations up to 3500/4000 m (Lüthi et al., 2019).
A benefit of the proposed downscaling approach is that is based on the truly local features, which were derived from 20 years of observations. In contrast to the final imputation step of SNCp50, which is based on a simple elevational dependence of snow cover, which can, in theory, directly be estimated from a low-resolution RCM signal. At least, the initial derivation of SNCp50 and the first imputation step can be assumed to provide downscaled estimates based on local features, while the 385 results of pixels that were subject to the final elevation imputation are more generalized. On the other hand, the pixel-by-https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License. pixel approach introduced artefacts at the low-resolution pixel boundaries (e.g., Fig. 5c). For the future, other downscaling techniques could be explored, such as analogue based or spatially explicit downscaling approaches based on principial components analysis.
One assumption in the downscaling is that the remotely sensed observations from MODIS are true, but these also have errors 390 and noise. Generally, accuracies in determining binary snow information (snow or land) are largely above 90% for MODIS (Parajka and Blöschl, 2006;Gafurov and Bárdossy, 2009). However, considerable uncertainty and lower accuracies were found for forested areas and locations affected by terrain shading (Notarnicola et al., 2013b). Specific to this study is the use of a cloud filtered product, which provides gap-free spatiotemporal series. The used filtering techniques resulted in only slightly lower overall accuracies of 91.5% compared to 93% for the original images (Matiu et al., 2020a). However, the 395 spatial and temporal filters that were applied to remove clouds might miss short snow episodes in low elevations and are difficult to validate in higher elevations, because of low ground station coverage. A pixel with erroneous information from MODIS will translate to an erroneous downscaled pixel, so relying on single pixels without consulting the spatial surrounding is not advised.
In addition, the downscaling assumes no land cover change, which might be problematic, for example, where the tree line 400 increases, and forests migrate to higher elevation. This comes on top to the already challenging estimation of snow cover fraction from remote sensing for forested areas. Under a warming climate, complex vegetation-snow interactions can occur, such as opposing effects on the interception and subsequent melting of snow in forests (DeBeer et al., 2021).

Comparison of downscaling to a dedicated snow model 405
For the Ötztal Alps region in Austria (Fig. 1), we compared results from bias correction and downscaling of RCM snow cover fraction (SNC) to running a dedicated snow model (AMUNDSEN), which has been forced by meteorology from RCMs. For the past period , the downscaling resulted in lower SNC than AMUNDSEN up to approximately 2000 m, similar SNC from 2100 to 2600m, and higher SNC for elevations above 2700 m ( Fig. 6 and S11). However, elevations above 2700 m are challenging to compare, since many pixels were removed from bias correction and downscaling 410 at these elevations because of snow accumulation issues and glaciers, while AMUNDSEN resolved the whole domain and explicitly considered ice-snow transitions. Consequently, comparisons above 2700 m are not based on the same pixels. such that decreases in SNC became stronger with increasing elevation. This elevation gradient in the RCM signal was also 420 present at the whole study region scale (Fig. S10b), but there downscaling showed consistently stronger decreases than bias correction only.
We assume this different behavior across elevations might be partly caused by elevation dependent warming, which implies stronger warming rates with increasing elevations, and which was observed in GCMs and RCMs at varying resolutions (Palazzi et al., 2019;Kotlarski et al., 2015;Lüthi et al., 2019). This might explain the elevation gradient in SNC changes 425 based on RCM data. While the projections from AMUNDSEN were also based on RCM meteorology, they underwent significant post-processing, including bias correction, spatial interpolation, and temporal disaggregation. Given this study's setup, it's not possible to disentangle how climate change signals and uncertainties flow through the modelling chain of both approaches with their different statistical post-processing and physical models. But we propose that such an assessment would be beneficial for highlighting important aspects of the modelling uncertainty of future mountain snow cover. A related 430 issue is that for bias correction of SNC, RCMs caused more of the overall variability than their driving GCMs, while in Hanzer et al. (2018), it was the opposite.
A further cause of the strong differences in SNC changes especially at lower elevations might be due to the consideration of forest snow processes in the AMUNDSEN simulations, where a canopy submodule accounts for the interception of snow by the trees -from where the snow can subsequently sublimate or melt without reaching the ground -as well as the 435 modification of the meteorological variables for sub-canopy conditions; for details see (Strasser et al., 2011). As the AMUNDSEN SNC results considered in this study only correspond to snow on the ground, this can cause differences with the RCM-based SNC changes, considering the large proportion of forested areas in the affected elevation bands (61 % forest coverage for elevations < 2000 m compared to only 2 % for elevations >= 2000 m). from RCMs. Shaded grey area (above 2700 m) indicates elevations, where >20 % of the pixels entering the average per elevation band were removed from MODIS and Downscaling but remained included in AMUNDSEN: these consist of glacierized pixels or pixels subject to snow accumulation in RCMs, while AMUNDSEN resolved the whole domain.

450
To conclude the comparison of bias correction and downscaling to using a dedicated snow model, Table 1 offers an overview of their main features. Both approaches enable to assess climate model uncertainty by using model ensembles, but both suffer from the potential need to extrapolate the RCM signal (surface meteorology or snow cover) beyond its elevation coverage. The main differences between the two approaches are in their spatial extent, spatial detail, and the representation of snow and ice processes. 455 Table 1. Non-exhaustive comparison of benefits (denoted with a +) and drawbacks (denoted with a -) of the two methods considered in this study. Abbreviations: regional climate model (RCM), snow cover fraction (SNC).  (Table S5). Their estimates were based on the Alpine3D snow model for 465 subregions of Switzerland, forced by RCM meteorology from the ENSEMBLES project, under the A2 emission scenario, which has less GHG concentrations at the end of the century compared to the RCP8.5 in this study. The differences in the reference values by elevation likely originate from the different study extent, while the differences in change estimates might partly be caused by different reference periods for the past. their domain includes much higher elevations due to the lower horizontal spacing, and the higher elevations showed stronger reductions in snow cover. 475

Conclusion
Bias correction of snow cover fraction from RCMs using aggregated MODIS remote sensing observations offers a promising approach to evaluate future changes of snow cover fraction in mountain areas under a climate model ensemble view. While limited by the resolution of RCMs, it offers consistent large-scale patterns of snow cover fraction for the past and future, and it is potentially applicable on a global scale. Consequently, it might be a viable alternative in remote or less monitored areas. 480 While snow cover fraction is probably not a priority for climate modelling groups to be made available, the proposed bias correction could benefit from an increasing ensemble of climate models being present. Regarding bias corrections methods, trend-preserving approaches, such as delta change or quantile delta mapping, were found superior to quantile mapping for snow cover fraction.
For the study region, which is approximately the Greater Alpine Region, results showed an overall reduction in snow cover 485 fraction for 2071-2100 compared to 2000-2020 of 17 % for RCP2.6 and 49 % under RCP8.5. However, strong elevational and seasonal dependencies of changes were found (Table S3 and S4, Fig. S5), with stronger decreases at middle elevations in winter (500-2000 m) and higher elevations in spring (above 1500 m). In addition, spatial patterns of change emerged, with stronger relative decreases in the south and west compared to north and east (Fig. 5), which are consistent to past trends of station observations of snow depth (Matiu et al., 2021). 490 The downscaling of RCM SNC with high-resolution MODIS observations falls under an "experimental" label. It suffers from many inadequacies, such as snow accumulation in RCMs, noise in observations, and is likely inappropriate for glacierized areas. However, it can provide auxiliary and high spatial resolution information while accounting for climate model uncertainty. The discrepancies to results from a dedicated snow model require further research before a final recommendation can be given. Potential usages of the downscaled information include hydrological studies or glacier 495 modelling studies that require snow line information. They might help in determining winter sport reliability, even though future assessments that do not account for technical snow are most likely not very useful (Spandre et al., 2019;Morin et al., 2021). Finally, downscaling approaches should be kept in mind considering the new generation of soon to be available highresolution (at or below 2 km) RCMs, for example, from CORDEX flagship pilot studies, together with long-term remote sensing observations at tens of meters scale, such as harmonized Landsat Sentinel series. 500 505 https://doi.org/10.5194/hess-2021-449 Preprint. Discussion started: 23 September 2021 c Author(s) 2021. CC BY 4.0 License.