Assessment of snow model uncertainty in relation to the effect of a 1 °C warming using the snow modelling framework openAMUNDSEN

Novel climate model data at the kilometer-scale, innovative downscaling techniques, sophisticated snow modelling frameworks, and increasing computational capacities are among the elements currently paving the way for a new phase in high resolution and physically based climate impact studies on snow hydrology in complex mountain terrain. While the assessment of climate model uncertainty is well established, the uncertainty arising from snow model selection typically receives far less attention. To investigate the uncertainty induced by the selection of the snow model configuration, we simulate the seasonal snow cover in the complex mountain area of the Berchtesgaden National Park mountains (Germany) under historical conditions (October 2013–September 2023) and for a 10 year period characterized by a 1 °C warming, using a large number of openAMUNDSEN snow model configurations (n=108) with degree-day as well as physically based snowmelt methods and varying land cover maps and spatial resolutions. The analysis of the resulting snow cover durations and snow disappearance days indicates that differences due to the choice of snowmelt method, land cover map and spatial resolution can be comparable in magnitude to the effect of a 1 °C warming, with uncertainties particularly pronounced in the forested areas and in the high elevations of the study area. Our results support the identification of critical snow model settings that need careful consideration, especially when employing energy balance instead of degree-day snow models to investigate climate change impacts on snow hydrology in complex mountain terrain.

1 Introduction

Climatic changes are expected to strongly affect the build-up and melt of seasonal snow covers in mountain regions. Because snow cover regulates the seasonal rhythm of mountain ecosystems and represents a key component of the hydrological cycle, a detailed knowledge of its past, current and future state is of great significance (Gobiet et al., 2014; Beniston et al., 2018; Matiu et al., 2021; Rumpf et al., 2022; Kotlarski et al., 2023; Vormoor et al., 2025). In addition to in-situ observations and remote sensing data, snow models are often used to assess the current state of the snow cover and the effects of climatic changes. Such impacts of climatic changes are typically evaluated by forcing the snow-hydrological models with bias adjusted scenario data from ensembles of climate models (e.g. Bürger et al., 2011; Marke et al., 2015; Vormoor et al., 2015; Muelchi et al., 2021). Thereby, empirical snow models using degree day factors (DDFs) are often used in climate change impact studies because they require little input data and are computationally very efficient. These temperature index (T-Index) approaches are typically optimized using seasonally and/or spatially varying DDFs (Hock, 2003; Zhou et al., 2021). Hanus et al. (2021), for example, use a semi-distributed hydrological model with an improved degree-day method including a seasonally variable DDF (Girons Lopez et al., 2020) to assess runoff changes in Alpine catchments in Austria. In the SNOWGRID-CL model, which was applied to the Ötztal Alps (Austria) presented in the climate change impact study by Kotlarski et al. (2023), an extended degree-day approach that accounts for air temperature and shortwave radiation was used (Olefs et al., 2020). The degree-day approach is also utilized in the semi-distributed HBV model (Bergström, 1995; Seibert and Bergström, 2022) or in the GR models (Coron et al., 2017) (e.g. Vormoor et al., 2015; Adhikari et al., 2022; Dakhlaoui et al., 2022; Ten Berge et al., 2025). To better describe rain-on-snow events, the mesoscale Hydrological model (mHM), which has been applied to assess changes in low flows and river floods at different warming levels (Thober et al., 2018; Marx et al., 2018; Rottler et al., 2021), accounts for energy input from liquid rain; DDFs are thereby increased depending on the amount of liquid precipitation (Samaniego et al., 2010).

In addition to degree-day approaches, energy balance based snow models, which explicitly represent the relevant energy and mass fluxes, are used to simulate the seasonal snow cover. These physically based modelling approaches require more input data and go along with a higher computational effort. Despite the increased complexity and effort, first studies have already employed energy balance based melt process description methods to assess the climate change impact on the seasonal snow cover. Hanzer et al. (2018), for example, use an energy balance based snowmelt modelling approach to assess cryospheric and hydrological impacts in the Ötztal Alps (Austria). They process an ensemble of climate model data using statistical bias correction methods as well as temporal disaggregation tools to conduct their energy balance based simulations. Marty et al. (2017) force the Alpine3D snow pack model with an ensemble of downscaled climate data for three emission scenarios in two Swiss catchments. A semi-distributed approach using representative elementary areas and the Crocus snow model (Brun et al., 2012; Vionnet et al., 2012) has been used to simulate the snow cover dynamics in French mountain massifs for historical conditions (Vernay et al., 2022) and to assess potential climate change impacts (Verfaillie et al., 2018; Kotlarski et al., 2023). A multi-physics ensemble approach for Crocus, which supports the representation of modelling errors, is presented in Lafaysse et al. (2017).

When designing a climate change impact study for complex mountain terrain, not only the selection of the snowmelt modelling method needs to be considered, also the need for lateral snow redistribution (Freudiger et al., 2017; Hanzer et al., 2016; Mott et al., 2018; Quéno et al., 2024) and snow-canopy processes (Pomeroy et al., 2012; Sun et al., 2022; Sanmiguel-Vallelado et al., 2022; Essery et al., 2025) have to be addressed. The selection of the model configuration should be tailored to the specific research goal. In practice, however, the data availability and computational resources are often limited. Currently, newly available ensembles of kilometer-scale climate data (e.g. Ban et al., 2021; Pichelli et al., 2021; Collier et al., 2024) are paving the way for a new generation of climate change impact studies in mountain areas. However, while the uncertainty in the climate scenario data is typically addressed using an ensemble of model runs, the snow model uncertainty receives only little attention. Further studies that quantify the effect of different snow model configurations are therefore required to support the development of robust modelling approaches for climate change impact studies in seasonally snow covered mountain areas.

In this study, we use the open source, intermediate complexity mountain snow model openAMUNDSEN (Strasser et al., 2024) to simulate the build-up and melt of the seasonal snow cover in the complex mountain terrain of the Berchtesgaden National Park (BGNP) under historical conditions (October 2013–September 2023) and for a 10 year period characterized by a 1 °C warming. The warming period is obtained using a stochastic block bootstrap resampler (climate generator), which is available as external pre-processing routine of openAMUNDSEN. We conduct snow simulations under historic and warming conditions using degree-day as well as energy balance based snow simulation routines, different land cover maps and grid geometries for the modelling with varying spatial resolutions (Fig. 1). For the historical period, model runs are evaluated against snow coverage maps derived from Sentinel-2 and MODIS data. For each model run, we calculate snow cover durations (SCDs) and snow disappearance days (SDDs) across the entire elevation range of the BGNP. Our main goal is to assess differences in the snow modelling results induced by different snow model configurations relative to the effect of a 1 °C warming. The results support the need for careful consideration in the choice of snow model type, input data, and spatial resolution. The simulated effects of a climate change signal need to be interpreted having these uncertainties in mind.

Figure 1 Schematic overview of the modelling set-up illustrating the snow cover simulations conducted with the mountain snow model openAMUNDSEN for the BGNP area under historical and 1 °C warming conditions using different snowmelt simulation methods, land cover maps and spatial resolutions.

2 Study area and data

We identified the BGNP as an ideal study site for our snow modelling experiment. The BGNP covers an area of 208 km² and is located in the state of Bavaria (Germany) in the North-East of the European Alps. It is the only German national park in the Alps and is characterized by a complex mountain landscape with an extreme topography covering elevations between 603–2713 m a.s.l. (meter above sea level) (Figs. 2 and 3). The open areas with bare rock and Alpine grassland in the higher elevations transition into coniferous, mixed and deciduous forest toward the lower elevations. Annual precipitation ranges from 1500 mm in the lower elevations to roughly 2600 mm in the peak regions, and the build-up and melt of a seasonal snow cover dominates the hydrology in the area (e.g. Strasser, 2008; Warscher et al., 2013; Storebakken et al., 2025). In cooperation with the Bavarian Avalanche Warning Service, the administration of the BGNP operates a dense network of automatic weather stations (AWSs) in the area. A detailed description of the station network including station names, locations, measured variables, example time series, and a description of the pre-processing of the AWS data is presented in Storebakken et al. (2025). The BGNP administration also provides a digital elevation model (DEM) in 10 m resolution, which we resample to resolutions between 20 and 1000 m (Fig. 3), and the BGNP habitat map (Fig. 2a). This habitat map was reclassified to seven classes, the three forest classes being the most important because forest types govern snow-canopy interactions and thus the entire seasonal dynamics of the snow cover in the forested areas. In addition to the BGNP land cover map, we process five other land cover maps for the area, all based on freely available and widely used datasets (Fig. 2b–f). Our aim is to capture the uncertainty that can be introduced into the distributed snow modelling with the selection of the land cover map. For comparability, we re-classify the data to similar land cover types. With regard to the forest representation, some data sources only have one forest class, whereas others distinguish different forest types and densities (Fig. 2). The land cover maps that we processed are based on the following data: (i) CORINE Land Cover 2018 obtained from the European Union's Copernicus Land Monitoring Service (https://doi.org/10.2909/960998c1-1870-4e82-8051-6485205ebbac, last access: 8 July 2025), (ii) OpenStreetMap (OSM) (https://www.openstreetmap.org/, last access: 8 July 2025), which we download using the R package “osmdata” (https://github.com/ropensci/osmdata, last access: 8 July 2025) (Padgham et al., 2017), (iii) WorldCover from the European Space Agency (ESA) (https://esa-worldcover.org/, last access; 8 July 2025) (Zanaga et al., 2022), (iv) Land Use/Cover Map for the European Alps developed by Marsoner et al. (2023) and (v) Land Cover Germany from the German Aerospace Agency (DLR) described in Weigand et al. (2020) (https://doi.org/10.15489/1ccmlap3mn39, last access: 7 July 2025).

Figure 2 Land cover maps for the BGNP based on the National Park habitat map (a; BGNP), CORINE land cover data (b; CORINE), OpenStreetMap data (c; OSM; https://www.openstreetmap.org/copyright, last access: 14 April 2026), ESA WorldCover data (d; WorldCover), the land use/land cover map for the European Alps developed by Marsoner et al. (2023) (e; EUSALP) and the Land Cover DE data from the German Aerospace Center (DLR) (f; LandCover-DE). The percentage coverage of each land cover class within 50 m elevation bands is shown right of each map.

Figure 3 Digital elevation models for the BGNP in different resolutions including the location of automatic weather stations (AWS): 20 m (a), 50 m (b), 100 m (c), 250 m (d), 500 m (e) and 1000 m (f). The elevation distribution within 50 m elevation bands is shown right of each map.

3 Methods

3.1 Intermediate complexity snow modelling

We use the open source snow-hydrological model openAMUNDSEN v1.0.4 to simulate the seasonal snow cover of the BGNP. openAMUNDSEN is a fully distributed snow model primarily designed to simulate the mass and energy balance of a seasonal snow cover in complex mountain terrain. On top of that, two different degree-day snowmelt modelling approaches are available as alternatives, namely the T-Index and the enhanced T-Index method. A comprehensive model description including these T-Index methods and all important process representations is presented in Strasser et al. (2024), the source code of the model is available on GitHub (https://github.com/openamundsen/openamundsen, last access: 7 July 2025). In the T-Index simulations presented here, we run the model in daily resolution with a constant DDF of 2.1 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. This DDF was selected based on comparisons between model results with satellite-based snow data (see Sect. 3.3). No automated calibration was conducted, but a robust value for DDF was selected following a series of model experiments. In addition to the T-Index approach, which only takes temperature into account, we also use the enhanced T-Index approach available within openAMUNDSEN, which also includes solar radiation (either measured radiation or potential clear-sky radiation) and modelled snow surface albedo (Pellicciotti et al., 2005). The incorporation of a radiation-driven melt component supports seasonally and spatially varying snowmelt dynamics and enables simulations to be run at the sub-daily scale. The enhanced T-Index snow model runs was conducted in a 3 hourly time step with a DDF of 1.3 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ and a short-wave radiation factor (RF) of 0.1 ${\mathrm{m}}^{\mathrm{2}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{W}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$. In openAMUNDSEN, the T-Index and enhanced T-Index snowmelt modelling approaches are combined with the “cryolayers” snow layering scheme, which enables to distinguish new snow, old snow, firn and ice. On the contrary, our 3 hourly openAMUNDSEN energy balance model runs engage a multilayer snowpack scheme following Essery (2015), distinguishing up to three snow layers (Strasser et al., 2024). For all model runs, we corrected precipitation for wind-induced measurement errors following the approach presented by Kochendorfer et al. (2017). The station-based meteorological forcing data was spatially interpolated using a combined lapse rate/inverse distance weighting scheme. First, lapse rates were applied to generate a distributed elevation-dependent field for each meteorological variable (i.e., the regression field). Next, residuals between the regression field and station observations were computed and interpolated to the grid using inverse distance weighting (IDW), yielding the residual field. The superposition of the regression and residual fields produces the final meteorological input grid. This interpolation procedure is performed at each model time step. We accounted for the lateral snow redistribution process with a simple snow redistribution parameterization using snow redistribution factors (SRF), calculated based on the concept of negative topographic openness (Yokoyama et al., 2002). An evaluation of this approach with airborne laser scanning data is presented in Hanzer et al. (2016) and a description of an elevation-dependent adjustment of the SRF, tailored for the BGNP and also used in this study, is available in Storebakken et al. (2025). We assessed SRFs for different resolutions and used them to redistribute the amount of solid precipitation from areas of erosion to areas of deposition (Fig. A1). All model input data and parameter settings used in this study are available in the BGNP model set-up data repository (https://doi.org/10.23728/B2SHARE.530A7560A736474599 69F5C21639E8CB, last access: 28 July 2025).

3.2 Stochastic climate generator

To assess snow model uncertainty in relation to a climate warming signal, we used the stochastic block bootstrap resampler (climate generator) from openAMUNDSEN to construct a 10 year period characterized by a 1 °C warming compared to the historical 10 year period. A detailed description of the climate generator is presented in Strasser (2008) and Strasser et al. (2024). The source code is available along with the openAMUNDSEN model (https://github.com/openamundsen/openamundsen-climategenerator, last access: 21 July 2025). The openAMUNDSEN climate generator is designed to construct time series with predefined mean temperature and precipitation. Therefore, the observational time series are split up into blocks of a defined length (one week in our case) and then – week by week – rearranged to obtain the new time series with the desired characteristics. For the selection of a week in the new time series, the algorithm cannot draw from all weeks of the year, but can only choose weeks observed up to two weeks before or after the week to be filled. The week 20 in the constructed time series, for example, can be selected from the weeks 18–22 in the observations. The width of this selection window is chosen by the user allowing a tailored compromise between a large data pool for selection and keeping seasonal characteristics consistent. Prior to the selection process, a random variability is added to the observations. Once an observed week has been selected by the algorithm to be used in the new time series, all measured variables (i.e., air temperature, precipitation, global radiation, relative humidity, and wind speed) are copied into the new data time series. This approach enables the generation of high temporal resolution (hourly in our case) and physically consistent station time series that directly can be used to force the snow model simulations. However, while stochastic block bootstrapping procedures can be very useful in the context of a snow model sensitivity study like the one presented here, we do not recommend their application in sophisticated climate change impact studies.

3.3 Model evaluation using satellite-based snow data

We evaluated the historical snow simulations using snow cover maps derived from Sentinel-2 and MODIS satellite data. The Sentinel-2 constellation allows for the determination of fractional snow cover (FSC) in high spatial resolution (i.e. 20 m×20 m). In this study, we use the Copernicus FSC obtained via the European Union's Copernicus Land Monitoring Service (https://doi.org/10.2909/3e2b4b7b-a460-41dd-a373-962d032795f3, last access: 30 July 2025). Copernicus FSC data is available from September 2016 onward. The available scenes for the study area were downloaded and pre-processed by cutting out the BGNP and re-sampling to the different model grids. We performed a pixel-by-pixel comparison for all scenes with less than 50 % cloud cover in the non-forested areas. The comparison of modelled and satellite-based data has four possible outcomes: true positive (TP), true negative (TN), false positive (FP) and false negative (FN). We determined the accuracy (ACC), overestimation error (OE) and underestimation error (UE) using the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\text{ACC}={\displaystyle \frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}}\text{(2)}& {\displaystyle }& {\displaystyle }\text{OE}={\displaystyle \frac{\text{FP}}{\text{FN}+\text{FP}}}\text{(3)}& {\displaystyle }& {\displaystyle }\text{UE}={\displaystyle \frac{\text{FN}}{\text{FP}+\text{FN}}}\end{array}$$

In comparison to, e.g., Hao et al. (2019) and Rottler et al. (2024), we adapted the OE and UE so that the errors are assessed as a fraction of the total amount of non-matching pixels. In this way, OE and UE sum up to 1 and deviations of the ACC from a perfect match can be attributed more directly to either an OE or UE. Alternative performance measures supporting snow model evaluation using satellite-based snow data, inter alia, are presented in Warscher et al. (2013) and Hofmeister et al. (2022). With regard to the model evaluation using FSC data, we present results on a monthly basis for February to July. Thereby we only take non-forested areas into account, as the satellite-based observation of snow on the ground inside forests is obstructed by the canopy (e.g. Xin et al., 2012; Bair et al., 2021; Keuris et al., 2023). The FSC and SWE values were translated to a binary classification using a threshold of 20 % for FSC and 20 mm for SWE. A sensitivity assessment with regard to this threshold selection is presented, e.g., in Rottler et al. (2024).

In addition to FSC maps, we use MODIS-based snow cover extent (SCE) to evaluate historical snow model simulations. We downloaded and processed the daily and cloud-free “Global SnowPack” product from the German Aeorospace Center (DLR) using DLR's GeoService (https://geoservice.dlr.de/, last access: 30 July 2025). This SCE product has a spatial resolution of 500 m. A description of this SCE product and examples of its applications in research studies are available in Dietz et al. (2012, 2013, 2014, 2015), Rößler et al. (2021) and Roessler and Dietz (2023). We calculated the seasonal course of the snow coverage in the non-forested areas of the BGNP for the 10 year period October 2013–September 2023 and compared the satellite-based estimates with openAMUNDSEN simulations.

4 Results

The evaluation of the openAMUNDSEN snow simulations in non-forested areas, using the satellite-based snow data, indicates that all model configurations can well capture the build-up and melt of the seasonal snow cover in the study area. A detailed description of the snow model evaluations is provided in the Appendix B. The pixel-by-pixel comparison of simulated snow with FSC from Sentinel-2 (Figs. B1 and B2) as well as the comparison with MODIS snow cover extent data (Figs. B3 and B4) confirm the model's ability to capture the seasonal snow cover dynamics in the BGNP.

The stochastic block bootstrap resampler (climate generator) generated a 10 year period characterized by 1 °C warmer temperatures (Fig. 4a). Due to the stochastic nature of the approach as well as the inherent natural variability of the observational data, not every month exhibits the exact same warming. Temperatures in December, for example, are on average 1.25 °C warmer, whereas March shows a more moderate warming of 0.61 °C. With regard to precipitation amounts, slight increases/decreases show up (Fig. 4). Again, theses differences reflect the stochastic nature of the method as well as the inherent variability in the data. We consider this variability in the signal appropriate for 10 year averages. In our experience, similar variabilities in the signal can show up in the comparison of 10 year periods of observations and data from climate models.

Figure 4 Monthly average temperatures (a) and precipitation totals (b) for the BGNP based on historical observations (October 2013–September 2023) and a 10 year period characterized by a 1 °C warming generated using the openAMUNDSEN climate generator. The average seasonal cycle of simulated snow water equivalent (c) and snow cover extent (d) as well as the average snow cover duration (e) and snow disappearance day (f) for 200 m elevation bands is assessed for all openAMUNDSEN model configurations (i.e. 108 for historic and for 1 °C warming conditions).

With regard to SWE, the snow model spread caused by the usage of different model configurations increases with progressing snow season and reaches its maximum between February and April. The average SWE in the BGNP can differ by more than 100 mm. Although the ensemble mean of the average SWE in warming conditions is below the ensemble mean of the historical simulations, there is a strong overlap of the shaded areas that indicated the spread caused by the usage of different snow model configurations (Fig. 4c). All model simulations under historical as well as warming conditions suggest that the BGNP is almost fully snow-covered in January and February (Fig. 4d). Furthermore, our snow simulations hint at an earlier disappearance of snow in the BGNP under a 1 °C warming by approx. one week. However, the strong overlap of the shaded areas indicates that the model spread in the second half of the main snowmelt season (May–June) is large compared to the magnitude of this shift. The model spread in this part of the year is characterized by values of approx. 50 mm for SWE and 12.5 % for the snow cover extent (Fig. 4c and d).

SCDs reach the highest values of approx. 200 d in elevation zones between 1800 and 2400 m (Fig. 4e). Our simulations suggest that a 1 °C warming can cause an average decrease in SCD of 20.5 d in the BGNP. For SDD, our modelling experiment suggests that changes due to a 1 °C warming can be small compared to the differences arising from the different model configurations (Fig. 4f). In general, our model results suggest that for SCD and SDD small differences between model configurations show up in elevations between 1600–2200 m. A larger spread in the snow model results can be detected below 1600 m and above 2200 m (Fig. 4e and f).

Our analysis of the effect of the spatial resolution suggests that results are largely consistent for resolutions between 20 and 250 m and up to elevations of 2000 m (Fig. 5). Simulation results for coarser spatial resolutions (i.e., 500 or 1000 m) can deviate from the ones obtained for higher resolution model runs. The median SCD for the elevation band 2200–2400 m simulated using the T-Index approach (Fig. 5a), for example, is 229, 218, 207 and 192 d using a spatial resolution of 20, 50, 100 and 250 m, respectively.

Figure 5 Simulated snow cover durations and snow disappearance days in the BGNP for different grid resolutions using the land cover map derived from the BGNP habitat data and the T-Index (a, b), enhanced T-Index (c, d) and energy balance (e, f) approaches available in the openAMUNDSEN snow model. Results are aggregated by 200 m elevation bands. Boxplots are displayed without whiskers and outliers.

The analysis of the model simulations with regard to the usage of different land cover maps indicates that the T-Index approach is largely insensitive to land cover (Fig. 6a and b). However, when using the enhanced T-index or energy balance based snowmelt approach, differences in SCD and SDD can show up in elevation zones with forest classes (Fig. 6c–f). The median SCD or SDD for elevations bands can be shifted by up to one week when using different land cover maps. In high elevations with no or only very small forested areas, the selection of the land cover map has no or only a negligible effect on the modelling results.

Figure 6 Simulated snow cover durations and snow disappearance days in the BGNP for different land cover maps using a 50 m spatial grid resolution and the T-Index (a, b), enhanced T-Index (c, d) and energy balance (e, f) approach available in the openAMUNDSEN snow model. Results are aggregated by 200 m elevation bands. Boxplots are displayed without whiskers and outliers.

Our model experiments indicate that a 1 °C warming can cause a reduction in SCD of more than 20 d in the BGNP (Fig. 7a, c and f). With regard to SDD, the change signal in our model experiment is smaller with values of a few days up to ten days (Fig. 9b, d and f). The changes induced by higher temperatures and assessed for elevation bands seem consistent for spatial resolutions between 20–250 m. Deviations between model configuration using different spatial resolutions mainly show up at coarser modelling resolutions (i.e., 500 or 1000 m) and for elevations below 800 m and above 2200 m. Our results also indicate that the magnitude of the difference can vary strongly within an elevation band. The length of the boxes, which cover 50 % of the cells within the elevation band, can easily span across a range of 10 d (Fig. 7).

Figure 7 Changes in simulated snow cover durations and snow disappearance days in the BGNP induced by a 1 °C warming using the land cover map based on the BGNP habitat map for different grid resolutions using the T-Index (a, b), enhanced T-Index (c, d) and energy balance (e, f) approach available within the openAMUNDSEN snow model. Results are aggregated by 200 m elevation bands. Boxplots are displayed without whiskers and outliers.

Our analyses of the effect of using different land cover maps on a potential climate change signal indicate that when using the T-Index approach, only small differences show up (Fig. 8a and b). When using the enhanced T-Index or energy balance approach, the usage of different land cover maps characterized by different extents and distributions of forest type classes can be stronger and alter results for elevation bands with forest (Fig. 8c–f).

Figure 8 Changes in simulated snow cover durations and snow disappearance days in the BGNP induced by a 1 °C warming for different different land cover maps using a 50 m spatial grid resolution and the T-index (a, b), enhanced T-index (c, d) and energy balance (e, f) approach available within the openAMUNDSEN snow model. Results are aggregated by 200 m elevation bands. Boxplots are displayed without whiskers and outliers.

Next, we directly compare results for SCDs and SDDs using different snowmelt simulation methods and assess results with respect to differences stemming from different aspect, land cover map and a 1 °C warming (Fig. 9). The aspects for the different model resolutions are presented in Appendix A (Fig. A2). Our results point out that the selection of the snowmelt modelling method can strongly influence SCD and SDD. Differences between the snowmelt modelling methods are most prominent below approx. 1600 m with simulations using the T-Index approach resulting in considerable shorter SCDs (Fig. 9a, c and e) and earlier SDDs (Fig. 9b, d and f). Energy balance based simulations, however, produce longer SCDs and later SDDs for this elevation zone, while results from the enhanced T-Index approach tend to be between those of the other two approaches. This signal reverses at high elevations (>2000 m) and the T-Index approach suggests longer SCDs and later SDDs than the other two snowmelt modelling approaches. Our modelling results indicate that at elevations below approx. 1600 m and above approx. 2200 m the differences induced by different snowmelt modelling approaches are stronger than the signal of change induced by a 1 °C warming (Fig. 9e and f). For elevation bands between 1600 and 2200 m, the results for SCDs and SDDs are more consistent and the change induced by a 1 °C warming exceeds the spread caused by different snowmelt simulation methods. Also the usage of different land cover maps can induce changes in SCDs and SDDs that can be within the same range as the change induced by a 1 °C warming (Fig. 9c and d), however, differences are less pronounced than differences resulting from the different snowmelt modelling methods, and are limited to elevations with differences in the forest representation in the land cover maps. While the enhanced T-Index and the energy balance approach can capture the differences in SCDs and SDDs depending on the exposition within elevation bands, the T-Index approach cannot account for such variabilities (Fig. 9a and b). Also the difference between northern and southern slopes is of the same magnitude as the signal induced by a 1 °C warming (Fig. 9).

Figure 9 Simulated snow cover durations and snow disappearance days in the BGNP using the T-index, enhanced T-index and energy balance approach available in the openAMUNDSEN snow model depending on aspect (a, b), land cover maps based on National Park habitat data (BGNP) and the land use/land cover map for the European Alps developed by Marsoner et al. (2023) (EUSALP) (c, d) and a warming of 1 °C (e, f). Results are aggregated by 200 m elevation bands. Boxplots are displayed without whiskers and outliers.

5 Discussion

The open source snow model openAMUNDSEN proved to be a versatile snow modelling framework offering a wide range of configuration options including degree-day and physically based snowmelt modelling approaches. The stochastic block bootstrap resampler (climate generator), an external pre-processing routine, is a useful extension of openAMUNDSEN and enables the generation of physically consistent meteorological times series with pre-defined warming trends. The incorporation of the climate generator into the modelling experiment enables us to analyse the snow model uncertainty in the context of a potential climate warming. While using the climate generator for our particular purpose – comparison of the results of different snow model configuration simulations in relation to a potential warming signal – is appropriate, we do not recommend employing a stochastic block bootstrapping procedure to produce climate change scenario data for impact studies, in particular, when the goal is the investigation of snow-hydrological extremes, as no new extreme events can be produced by simple re-organization of observational time slices. Moreover, it needs to be noted that the climate generator cannot adequately produce elevation- and season-dependent features of long-term climatic changes.

In our study, we compare averages of meteorological and snow-related variables for two 10 year periods. The average monthly temperatures and precipitation totals as well as the average seasonal cycles indicate that the inter-annual variability of the data is affecting the averages investigated (Fig. 4a–d). Our analysis highlights the need for caution when quantifying climate change signals obtained using such short time frames. However, in our study, we deliberately decided to compare 10 year periods, as this is the type of data currently being available from ensemble climate simulations at kilometer-scale resolution. These high resolution climate model runs better capture the topographic forcing of precipitation and allow for an explicit, physical description of deep convection without having to use parameterization schemes (Lucas-Picher et al., 2021). Ban et al. (2021), for example, conduct the first multi-model ensemble climate simulations for the greater Alpine region at the kilometer-scale for the ten year period 2000–2009. Multi-model 10 year kilometer-scale scenarios for the greater Alpine domain are presented in Pichelli et al. (2021). First case study simulations and results for a hydrological year for the “Third Pole”, using a multi-model and multi-physics ensemble at the kilometer-scale, are presented in Prein et al. (2023) and Collier et al. (2024), respectively. On the contrary, high resolution regional climate model runs covering longer time frames (e.g. 30 years or more) often are limited to a single model only (e.g. Warscher et al., 2019; Collier and Mölg, 2020). Furthermore, depending on the complexity of the terrain, resolutions below 1 km are required to adequately represent the processes that cause the spatial variability of the simulated snow cover, such as lateral redistribution of snow by wind or snow-canopy interaction. Consequently, further dynamical or statistical downscaling of kilometer-scale climate model data is required to assess changes in the seasonal snow cover dynamics in complex mountain terrain such as the BGNP. Recent studies have demonstrated that the combination of high resolution weather model output, e.g. from the Weather Research and Forecasting (WRF) model, with physically-based snow energy balance models is feasible and can provide new insights into snow dynamics at the regional scale (e.g. Corripio and López-Moreno, 2017; Liu et al., 2019; Alonso-González et al., 2021; Raparelli et al., 2023; Zhou et al., 2023). Ongoing improvements of the snow process descriptions in the weather models may eventually make the coupling to an additional snow model unnecessary in the future (e.g. Luo et al., 2021; Saigger et al., 2024).

The recent progress in climate and snow modelling research suggests that, step by step, computational hurdles will be overcome, and that future climate change impact studies using models that better capture the topographic forcing of precipitation and explicit treatment of deep convection, and physically based snow models including lateral redistribution and snow canopy processes can adequately describe the snow cover dynamics in complex mountain terrain for potential future warming conditions. In this study, we show that the selection of the snowmelt simulation method is a very important aspect to consider when conceptualizing climate change impact studies in complex mountain areas. Our modelling experiments indicate that the snow model uncertainty can be of the same magnitude or even larger than changes in the snow cover induced by a 1 °C warming (Fig. 9e and f). A particular focus is required for forested areas, as additional snow-canopy processes need to be taken into account. In contrast to degree-day approaches, which usually only consider temperature and a calibrated degree-day factor and are largely insensitive to different land cover types, energy balance snow cover models do react sensitively to the choice of the land cover map and the forest representation therein (Fig. 6), given the meteorological forcing is properly adjusted to inside-canopy conditions. The comparison of model simulations with different spatial resolutions suggests that for a robust quantification of changes in snow patterns in the highest elevations, spatial resolutions considerably below 1 km are required (Fig. 7). In our simulations, we use a simple snow redistribution parameterization to account for the lateral snow redistribution processes (Strasser et al., 2024) (Fig. A1). More process-oriented approaches incorporating wind fields and explicit representation of gravitational snow redistribution processes exist and need to be considered in future developments of the openAMUNDSEN model (e.g. Warscher et al., 2013; Freudiger et al., 2017; Quéno et al., 2024).

Our study contributes to the broader research objective of quantifying snow model uncertainty and follows up on the investigations by Essery and Etchevers (2004), Essery et al. (2013), and Günther et al. (2019, 2020) who assessed uncertainties in snow model results arising from snow model complexity, parameter errors and errors in the forcing data. The effect of process representations on the results are also presented in Magnusson et al. (2015). Particular consideration of snow in forested area was given in the second phase of the Snow Model Intercomparison Project (SnowMIP2) where output from 33 different models was compared (Essery et al., 2009). These studies highlight that the accuracy of the forcing data is more critical than parameter selection, while the type of the model seems of less importance. In addition, human errors are an important aspect affecting snow model performance (Menard et al., 2021). Our snow modelling experiments complement existing research by defining a large set of configurations for distributed snow model simulations in complex mountain terrain with open and forested areas, resulting in a large number of openAMUNDSEN model configurations (n=108) with both degree-day and physically based snowmelt simulations, varying land cover maps and different spatial resolutions. The assessment of results under a 1 °C warming allows for a quantification of snow model uncertainty relative to a potential warming signal.

As many seasonally snow covered areas are also forest covered, further research is needed to better understand the effect of the snow-canopy processes in the context of a changing climate. A direct coupling of forest and snow models will thereby be an important next step in this regard, since not only climate will change, but also species composition and hence the snow hydrological processes inside the canopy. The complex mountain research site BGNP seems to be an ideal study site for such a task, given the dense measurement network and the existing expertise in forest processes modelling (e.g. Albrich et al., 2022; Senf et al., 2017; Thom et al., 2022; Dollinger et al., 2023; Braziunas et al., 2024) as well as snow modelling (e.g. Strasser, 2008; Strasser et al., 2011; Bernhardt et al., 2012; Kraller et al., 2012; Warscher et al., 2013; Storebakken et al., 2025). At the same time, however, the evaluation of snow models inside forests remains challenging, particularly, as forest canopies obstruct the satellite-based detection of snow on the ground. An alternative approach, the evaluation of openAMUNDSEN snow modelling results inside the forest using a network of microclimate sensors, is presented in Storebakken et al. (2025). Further research will support the development of new approaches for the evaluation of distributed snow models applied to forested areas, e.g. with respect to plant physiological simulation of LAI or improved scaling function for the meteorological variables from the location of the measurement to inside-canopy conditions. In our view, such new techniques for snow model evaluation are required to prepare a new generation of distributed and physically-based snow models for the application in the climate change impact context. There is no physical bases for assuming that DDFs used to simulate historical conditions are valid under future climatic conditions, and recent investigations of DDFs based on energy flux components (Ismail et al., 2023) show that DDFs cannot be treated as constant factors. They are changing with progressing melt season and need to be specifically adapted. Also our results suggest that degree-day snowmelt modelling approaches have their limitations, i.e. the DDF needs careful calibration and varies with elevation. We found that in higher elevations the best DDF is larger than in the valley region, where intermittent melting of the shallow snow cover in early winter can be very fast. Accordingly, the most consistent results and only small differences between model configurations show up in elevations between 1600–2200 m (Figs. 4e, f and 9e, f). The spatially and seasonally variable snowmelt dynamics from the enhanced T-Index approach seem to improve results by accounting for aspect- and season-dependent variations in radiation. However, the differences in SCDs and SDDs between the enhanced T-Index and the energy balance simulation runs can be substantial (Fig. 9). At lower elevations, the different forest representation in each particular land cover map causes additional uncertainties in the results, while SCDs and SDDs in the highest elevations react very sensitively to the selection of the spatial resolution (Figs. 6 and 5). Results are only consistent for spatial resolutions between 20 and 250 m and at elevations up to 2000 m.

6 Conclusions

openAMUNDSEN is a versatile snow modelling framework that supports distributed snow cover simulations in complex mountain areas, including snow-canopy interaction and lateral snow redistribution, using degree-day as well as energy balance snowmelt methods. In this study, we employ a large number of openAMUNDSEN model configurations (n=108) characterized by different snowmelt modelling methods, land cover maps and spatial resolutions to simulate the seasonal snow cover in the complex mountain terrain of the BGNP for the historical period October 2013–September 2023 and a 10 year period characterized by a 1 °C warming. The data describing the warming period was obtained using a stochastic block bootstrap resampler (climate generator), which is available as external pre-processing routine of openAMUNDSEN. The incorporation of the climate generator into our modelling experiment enabled us to assess snow model uncertainty in relation to the effect of a potential warming signal. The complex and diverse mountain landscape of the BGNP proved to be an ideal study area for snow modelling research and the presented investigation of snow model uncertainty. Our analyses suggest that differences in the resulting SCDs and SDDs, arising from the selected snowmelt simulation method, the used land cover map and the chosen spatial resolution can be comparable to the impact of a 1 °C warming. The uncertainties in the results are pronounced in the forest covered areas and in the high elevations of the study area. Our results emphasize the importance of evaluating the uncertainty of a given snow modelling setup when quantifying the effect of a warmer climate. Sources of uncertainty like the selection of the snowmelt method, land cover map, or the chosen resolution can have effects on the results in the same order of magnitude as the climate change signal. Future snow modelling studies investigating the effect of climatic changes on the seasonal snow cover in complex mountain terrain can build on our findings by carefully selecting snow model configurations and related uncertainties.

Appendix A: Snow redistribution factors and aspects

Figure A1 Snow redistribution factors for the BGNP for different spatial resolutions: 20 m (a), 50 m (b), 100 m (c), 250 m (d), 500 m (e) and 1000 m (f). The frequency of different snow redistribution factors is shown in the panel to the right of each map.

Figure A2 Aspects for the BGNP for different spatial resolutions: 20 m (a), 50 m (b), 100 m (c), 250 m (d), 500 m (e) and 1000 m (f). The frequency of different aspects is shown in the panel to the right of each map.

Appendix B: Model evaluation using satellite-based snow data

The pixel-by-pixel comparison based on FSC maps derived from Sentinel-2 images indicates an average accuracy for the months February to July between 0.82 and 0.93 (Fig. B2). Our results indicate that the coarser the spatial resolution of the model run and the re-sampled satellite image, the higher the accuracy obtained by this pixel-by-pixel comparison (Figs. B1 and B2). Using the energy balance approach, an overestimation of the snow cover shows up at the beginning of the snowmelt season. The deviations between simulated and observed snow cover at the end of the snowmelt seasons, however, are due to an underestimation error, i.e., too little snow covered area in the model simulations (Fig. B2c1–c6). The comparison of the simulated build-up and melt of the seasonal snow cover with MODIS snow cover extent data for the 10 year period October 2013–September 2023 indicates a good agreement between simulated and satellite-based data (Figs. B3 and B4). Looking at the evolution of the seasonal snow cover as the extent over the entire non-forested area, only very little differences show up between the different spatial resolutions, and the snow cover extent curves largely overlap for all snowmelt modelling approaches (Fig. B4). For individual years (e.g., 2015 and 2016), the empirical snowmelt modelling approaches seem to produce a too late melt of the snow cover in the high elevations (Fig. B4a and b). The calculated NSE values range from 0.86–0.93.

Figure B1 Simulated snow water equivalent (SWE) using the energy balance approach available in openAMUNDSEN and satellite-based fractional snow cover (FSC) maps derived from Sentinel-2 data for the BGNP on 23 April 2020 for different grid resolutions: 20 m (a), 50 m (b), 100 m (c), 250 m (d), 500 m (e) and 1000 m (f). The accuracy (ACC) assessed comparing the simulated SWE with the corresponding re-sampled satellite-based FSC observation is noted on the top-right of each SWE map.

Figure B2 Evaluation of openAMUNDSEN snow simulations in different spatial resolutions for non-forested areas in the BGNP between February and July using fractional snow cover (FSC) maps derived from Sentinel-2 data by means of accuracy (ACC), overestimation error (OE) and underestimation error (UE) for images with not more than 50 % cloud coverage. The average ACC for all scenes is noted on the top-right of the plot panels.

Figure B3 Simulated snow water equivalent (SWE) using the T-Index (a), enhanced T-index (b) and energy balance (c) snowmelt simulation approaches available in the openAMUNDSEN snow model with a spatial resolution of 100 m and satellite-based snow cover extend maps (SCE) derived from MODIS snow cover data (d) for the BGNP and the three selected days: (1) 16 April 2020, (2) 1 May 2020 and (3) 15 April 2020. The total snow covered extent (SCE) is noted on the top-right of each map.

Figure B4 Evaluation of the openAMUNDSEN snow simulations for non-forested areas of the BGNP conducted using the T-index (a), enhanced T-index (b) and energy balance (c) snowmelt simulation approaches available in the openAMUNDSEN snow model by comparing the evolution of the simulated snow cover extent with the snow cover extent from MODIS satellite data for the time period October 2013–September 2023. The Nash–Sutcliffe-Efficiency (NSE) for model runs with different spatial resolutions is noted on the top-right of each panel.