A high-resolution snow dataset for Switzerland (2016–2025) combining physics-based simulations and in situ observations

We present a high-resolution snow dataset that provides daily estimates of snow depth, snow water equivalent, snow cover fraction, and snowmelt runoff for Switzerland and hydrologically connected bordering regions, covering water years 2016–2025. The dataset is based on fully distributed simulations at 250 m resolution using the multi-layer, physics-based snow model FSM2OSHD, operated by the Swiss Operational Snow Hydrological Service. To capture the high spatial heterogeneity of snow cover dynamics in complex mountainous terrain, the modeling framework combines dedicated dynamical and statistical downscaling of numerical weather prediction data with the upscaling of hyper-resolution terrain, forest, and light-availability datasets, explicitly accounting for subgrid variability. The particle filter-based assimilation of in situ snow depth observations from 444 monitoring stations across the domain dynamically corrects spatiotemporal error patterns in the meteorological forcing data. This approach ensures consistent input data quality over the entire 10-year period and mitigates potential discontinuities caused by changes within the numerical weather prediction system. Example applications demonstrate the dataset's ability to capture regional and interannual variability of snow water resources, snow cover extent, and snow duration. With 10 years of physically consistent estimates at high spatial and temporal resolution, this dataset represents, to our knowledge, the most accurate and comprehensive record of snow cover dynamics for Switzerland to date. It expands the snow data record for the European Alps and bridges the gap between coarse global reanalyses and detailed local observations. The dataset is publicly and freely available providing a valuable resource for a wide range of scientific and applied studies in hydrology, ecology, climate, and cryospheric research.

1 Introduction

Seasonal snow drives numerous hydrological and ecological processes (Han et al., 2024; Hale et al., 2023; Slatyer et al., 2022) and affects many socioeconomic aspects (Rasul and Molden, 2019). Snow-cover extent, as well as timing and intensity of snowmelt, have a direct impact on avalanche and flood hazards (Musselman et al., 2018; Li et al., 2019; Eckert et al., 2024), freshwater availability (Musselman et al., 2021; Siirila-Woodburn et al., 2021), hydropower production (Magnusson et al., 2020; Hou et al., 2025), and winter tourism (Moreno-Gené et al., 2018; Morin et al., 2021). However, estimating snow water resources in mountainous regions is particularly challenging due to the substantial spatial variability of the terrain and the snowpack (Grünewald et al., 2010; López-Moreno et al., 2013; Mott et al., 2018), the lack of accurate distributed measurements (Dozier et al., 2016; Gascoin et al., 2024), and the uncertainties inherent in snowpack estimated from numerical models (Menard et al., 2021).

The mountain snowpack remains severely undersampled despite continuous monitoring efforts (Largeron et al., 2020). Detailed, high-resolution hydrometeorological and snow datasets exist for small catchments, such as the Dischma catchment in Switzerland (Magnusson et al., 2025), the Izas catchment in the Spanish Pyrenees (Revuelto et al., 2017), or the Johnston Draw catchment in Idaho, USA (Godsey et al., 2018), but their spatial coverage is limited. Similarly, airborne lidar- or photogrammetry-based snow depth maps are available only for a few select catchments and at specific times during the season (Painter et al., 2016; Bührle et al., 2023). For larger scales and with roughly weekly temporal resolution, spaceborne optical sensors provide observations of snow cover extent and snow cover fraction (SCF) (Gascoin et al., 2019), but cannot provide direct information on snow depth or snow water equivalent (SWE). The NorSWE dataset (Mortimer and Vionnet, 2025) compiles SWE observations from more than 10 000 locations across the Northern Hemisphere over three decades, with the Alps represented by manual point measurements from 11 sites in Switzerland only. Spatially explicit snow depth observations from Sentinel-1 retrievals are possible under dry-snow conditions (Lievens et al., 2019, 2022), but truly reliable, high-resolution, spatiotemporally continuous SWE estimates remain elusive (Gascoin et al., 2024; Mortimer et al., 2020). While numerical models can provide such continuous estimates at any desirable spatial and temporal resolution, they are subject to inherent uncertainties in the parametrization and forcing data (Günther et al., 2019; Menard et al., 2021).

Reanalysis products provide estimates of past states by constraining numerical simulations through the assimilation of observational datasets. While global reanalyses within numerical weather prediction (NWP) systems, such as ECMWF's ERA5 (Hersbach et al., 2020), are widely used for climate monitoring, their snow-related variables offer only low resolution and accuracy in mountainous regions (Li et al., 2022; Monteiro and Morin, 2023; Kouki et al., 2023). Fiddes et al. (2019) presented an efficient method that couples sub-grid clustering of complex terrain, downscaling of global meteorological reanalysis data, and the assimilation of spaceborn SCF observations, to enable high-resolution ensemble-based snow reanalyses in mountain regions. Other recent efforts in the snow modeling community have produced detailed and dedicated long-term snowpack reanalysis datasets: the modeling and data assimilation system SNODAS provides daily snowpack and precipitation data at 1 km resolution over the contiguous US from 2003 onward (Barrett, 2003); daily estimates of SWE and SCF based on the assimilation Landsat SCF observations into a land surface model, coupled with a snow depletion curve, are available for the Sierra Nevada (Margulis et al., 2016) and the western US (Fang et al., 2022), at a resolution of 90 and 500 m, respectively; a similar dataset based on the joint assimilation of Landsat and MODIS SCF observations is available for High Mountain Asia (Liu et al., 2021); a daily 10 km gridded snow depth and SWE estimates over the Iberian Peninsula since 1980 are available based on the physics-based snow model FSM forced with downscaled ERA-interim data (Alonso-González et al., 2018, and references therein); daily 1 km gridded estimates of snow depth and snow cover duration between 1961–2020 are available for Austria, based on simulations of the SNOWGRID model in a climate configuration forced with gridded meteorological observations (Olefs et al., 2020); the IT-SNOW reanalysis provides daily estimates of snow states for Italy from water year 2010 onward by combining model simulations with in-situ and spaceborne osbervations (Avanzi et al., 2023); the S2M meteorological and snow cover reanalysis combines the meteorological analysis SAFRAN and the high complexity snow model Crocus within the SURFEX/ISBA land surface model, covering the semi-distributed massifs of the French Alps, Pyrenees, and Corsica from 1958 onwards (Vernay et al., 2022, and references therein).

Here, we present a continuous 10-year snow reanalysis dataset for Switzerland and hydrologically connected bordering regions, produced within the near-real-time modeling framework of the Swiss Operational Snow Hydrological Service (OSHD). The dataset provides daily estimates of snow depth, SWE, SCF, and snowmelt runoff by combining high-resolution simulations from the intermediate-complexity, physics-based snow model FSM2OSHD (Mott et al., 2023) with in situ snow depth observations from 444 stations (Oberrauch et al., 2024). To capture the high spatial variability of snow processes in complex Alpine terrain, the modeling chain employs dedicated dynamical and statistical downscaling of NWP forcing data to a spatial resolution of 250 m, solves the coupled mass- and energy-balance equation for multiple numerical snow layers, and explicitely accounts for differences in the atmospheric and snowpack procsesses of open, forested, and glaciated areas (Mazzotti et al., 2021). The simulations are based on upscaled versions of the most accurate hyper-resolution datasets currently available, including a 10 m digital elevation model and terrain surface model (swissALTIRegio, 2025; swissSurface3D, 2025), as well as a 1 m canopy height model and light availability maps that resolve terrain and vegetation shading down to individual trees for every hour of the year (Webster et al., 2025). Temporal consistency across the whole period is ensured through an assimilation scheme that homogenizes meteorological inputs despite changes in the source and processing level of the forcing data (Oberrauch et al., 2025; MeteoSwiss, 2024). The model has been tuned and validated continuously over the past decade against snow depth, SWE, SCF, and new snow observations from a dense station network (Mott et al., 2023; Cluzet et al., 2024; Haagmans et al., 2026).

To our knowledge, this dataset represents the most accurate and comprehensive snow reanalysis for Switzerland to date. By providing 10 years of physically consistent estimates of snow cover dynamics at a subkilometric resolution across a large part of the European Alps, it contributes to the snow database of a region characterized by highly complex terrain and diverse hydroclimatic regimes. The dataset helps bridge existing scale and knowledge gaps between coarse global products (e.g. Hersbach et al., 2020) and detailed local observations (e.g. Magnusson et al., 2025) and provides a robust foundation for hydrological, ecological, and cryospheric analyses (e.g. Jenicek et al., 2016; Xie et al., 2018; Floriancic et al., 2020; Publications Office of the European Union, 2023). The data is publicly and freely available under https://doi.org/10.5281/zenodo.17313889 (Oberrauch et al., 2026), offering substantial value to a wide range of scientific and applied users.

2 Data and model chain

The basis for the presented dataset is the fully distributed, physics-based, multi-layer snow model FSM2OSHD (Mott et al., 2023), operated by the Swiss Operational Snow Hydrological Service (OSHD). Forcing data from the Federal Office of Meteorology and Climatology, MeteoSwiss, are debiased and downscaled to the 250 m model resolution and dynamically corrected through the assimilation of in situ snow depth observations (Oberrauch et al., 2024). The following section provides a brief overview of the study area, the datasets used, the modeling chain, and the particle filter (PF) based assimilation scheme.

2.1 Study area

The presented dataset spans Switzerland and hydrologically connected neighboring regions of Austria, France, Germany, Italy, and Liechtenstein (Fig. 1). The model domain has a latitudinal and longitudinal extent of 272 and 365 km, respectively. It comprises 928 155 grid cells at a spatial resolution of 250 m, for a total area of over 58 000 km². Model elevations range from 180 to 4750 m, with 50 % of the grid points located above 966 m. Roughly 40 % of the domain is forested, with forest cover extending up to 2400 m, although less than half of the forest is situated above 1000 m (Szerencsits, 2012; Haagmans et al., 2026).

Figure 1 Terrain representation of the model domain based on the 250 m resolution digital elevation model. Major water bodies are shown in grey, the Swiss national border is outlined in black, and a shaded relief is used as a background map. The bottom-right inset panel displays the areal distribution of the domain across 250 m elevation bands, with forest and non-forest fractions shown in green and beige, respectively. The top-right inset panel shows the 16 hydrological units, collectively referred to as the Swiss Alps, used for validation against spaceborne SCF observations.

Switzerland lies in Central Europe and spans both sides of the European Alps and the adjacent lowlands, each characterized by distinct hydroclimatic regimes. The Alpine Main Ridge forms a climatic divide between Mediterranean air masses to the south and Atlantic air masses to the north (Gubler et al., 2023). Precipitation is generally high in the Alps, the Alpine foothills, and the Jura mountains (in the northwest of the country), with annual totals of around 2000 mm. In contrast, certain inner-Alpine valleys are sheltered from moist air masses and are therefore comparatively dry, with annual totals of about 700 mm or less (MeteoSwiss, 2025). Figure 2 displays average monthly air temperature and solid precipitation aggregated into four elevation bands, based on the forcing data over the presented 10-year period (with the interannual variability indicated by shadings).

Figure 2 Monthly climatology of near-surface air temperature (left) and snowfall sums (right) stratified by elevation bands over the period between September 2015–August 2025, based on the meteorological forcing data. Results are averaged over all model grid points within the respective elevation bands. Solid lines represent mean values, while the shading indicates the standard deviation across the 10 seasons.

2.2 Used datasets

MeteoSwiss officially replaced its operational forecasting model in June 2024 (MeteoSwiss, 2024), shifting from COSMO (Baldauf et al., 2011) to ICON (Zängl et al., 2015) as the dynamical core. Moreover, ICON data was reprocessed for a transitional period between 2020–2024. ICON is the higher-resolution successor to COSMO, and provides an improved representation of terrain-induced variability (Zängl et al., 2015; MeteoSwiss, 2023). Hence, for the water years 2016–2020, we used analysis data from COSMO, while from 2021 onward, ICON data were employed at different processing levels: (re-)forecast data for the water years 2021–2023 and analysis data from 2024 onward. Oberrauch et al. (2025) demonstrated that correcting spatiotemporal error patterns in the forcing data through the assimilation of point observations homogenizes snow model performance. Thereby, potential temporal discontinuities arising from changes in the forcing data sources are mitigated, as detailed below in Sect. 2.4.

The hourly NWP data are downscaled from 1 km to the 250 m model resolution using various statistical and dynamic downscaling schemes (outlined below in Sect. 2.3), all based on a detailed land-use and land-cover datasets. The basis is the swissALTIRegio digital elevation model (DEM) with a resolution of 10 m (swissALTIRegio, 2025), used to compute the 250 m model DEM, as well as the topographic position index, slope, aspect, and subgrid terrain variability for each grid cell. The CORINE land cover dataset (European Union's Copernicus Land Monitoring Service information, 2020) was used to distinguish between open, forested, and glaciated areas, which are handled differently by the FSM2OSHD model (outlined below in Sect. 2.3). Large water bodies where snow does not accumulate are excluded from the domain, indicated as grey areas in Fig. 1 and subsequent figures.

The SwissRad10 dataset provides domain-wide estimates of direct and diffuse radiation at 10 m resolution across Switzerland, accounting for vegetation and terrain shadowing as well as the sky view fraction (Webster et al., 2025). The dataset is derived from high-resolution airborne lidar data that resolves individual trees (swissSurface3D, 2025), using the Canopy Radiation Model (Webster et al., 2023). It provides a terrain-only scenario, as well as leaf-on and leaf-off canopy scenarios (of which only the latter is used in FSM2OSHD) at hourly resolution over a full annual solar cycle.

The assimilation scheme (see Sect. 2.4) uses snow depth observations from a network of 444 monitoring stations across Switzerland and neighboring regions in Germany, Austria, Liechtenstein, Italy, and France, spanning an elevation range from 230 to 2950 m (with 50 % above 1300 m). All stations are in flat, open areas, while locations known for unrepresentative snow conditions (e.g. due to strong snow drift) are excluded. The data are manually quality-controlled, including the removal of obvious outlier observations and the filling of small data gaps (informed by observations from neighboring stations). Snow depth observations are reported daily, with less than 2 % of station-days missing over the 444 stations and 10 years.

Spaceborne SCF maps derived from Sentinel-2 data, provided by the AlpSnow science activity within ESA's Alpine Regional Initiative (European Space Agency, 2025), are used for an independent evaluation of the dataset. The SCF retrievals are based on multi-spectral unmixing of Level-1C reflectance data for improved snow detection in shaded areas (Keuris et al., 2023), while clouds, glaciers, forests, and urban areas are masked out to mitigate inherent limitations of optical sensors, following the methods of Cluzet et al. (2024). As these data are available to us only for the water years 2018–2025, they have not been assimilated in order to avoid temporal discontinuities in the final dataset (as discussed in Sect. 4).

2.3 The OSHD modeling chain

FSM2OSHD is a multi-layer, physics-based model solving the coupled mass and energy balance for individual numerical snow layers at an hourly resolution, without directly accounting for snow microstructures and metamorphism (Mott et al., 2023). FSM2OSHD was originally based on the Flexible Snow Model (FSM2, Essery, 2015; Essery et al., 2025), incorporating additional process-based refinements adapted and tuned for the application within the OSHD.

The implemented snowpack process parameterizations (FSM2OSHD code available at https://doi.org/10.5281/zenodo.20410047, Mott et al., 2026) are briefly outlined below; please refer to Mott et al. (2023, and references therein) for more detailed information. Fresh snow density in the model is estimated from air temperature and wind speed during snowfall (Vionnet et al., 2012), tuned against snow board measurements, which capture the depth and weight of snow accumulated over a period of 24 h (Fierz et al., 2009; Haberkorn, 2019). Since such measurements inherently include initial settling, the tuning accounts for both fresh snow density and compaction over the first 24 h. The settling of the snowpack and the associated changes in snow density are computed via a viscosity-based overburden scheme adapted from Vionnet et al. (2012). Increasing snow weight compresses the underlying layers depending on the layer's viscosity, which, in turn, depends on layer density, temperature, and liquid water content. Thermal conductivity between snow layers is diagnosed from snow density (Douville et al., 1995; Essery, 2015), while the ground heat flux is estimated using a five-layer soil model (Cox et al., 1999; Essery, 2015). Turbulent fluxes at the snow–atmosphere boundary are parameterized following Essery (2015, Sect. 2.3.4), with transfer coefficients adjusted according to the Richardson number (Ri, Louis, 1979) and capped at Ri=0.2 to prevent an unrealistic stability-induced shutdown of the turbulent exchange. Liquid water retention in each layer follows a bucket-storage approach, with the storage capacity depending on the snow density (Anderson, 1976; Boone and Etchevers, 2001). Broadband snow albedo is computed with a snow-age-dependent decay function, following a linear rate for cold snow and an exponential rate for melting snow (Douville et al., 1995). Albedo values are reset to their maximum fresh-snow values after a minimum of 10 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ new snow accumulation over 24 h. Additional tuning accounts for bare ground piercing through thin snowpacks, as well as aspect- and slope-dependent differences in the albedo decay based on a comparison with spaceborne snow wetness and SCF observations (Cluzet et al., 2024). SCF in open, non-forested areas is estimated by tracking the seasonal evolution of snow depth and SWE depending on subgrid terrain variability (Helbig et al., 2021), while a simpler hyperbolic tangent model is applied within forests (Essery, 2015).

FSM2OSHD differentiates between open, forested, and glaciated terrain by simulating the snow cover separately for each land cover type and aggregating the results as a weighted average according to the respective land cover fractions within each grid cell. Glaciated terrain is simulated analogously to open terrain but with a modified ground heat flux and roughness length. For forest-covered grid cells, FSM2OSHD explicitly accounts for key snow–canopy and canopy–atmosphere interactions, including snowfall interception, snow unloading and sublimation from the canopy, shortwave radiation transmission, longwave radiation enhancement, and wind attenuation (Mazzotti et al., 2020a, 2021; Essery et al., 2025). Forest processes and vegetation shading extend beyond the forest edges, thereby influencing adjacent open terrain.

The hourly NWP data are bias corrected and downscaled to the 250 m model grid and station locations (see Sect. 4 and Table 1 in Mott et al., 2023) by statistically downscaling the wind fields (Winstral et al., 2017), dynamically downscaling the radiation input (Jonas et al., 2020; Webster et al., 2025), and linearly interpolating air temperature, relative humidity, and precipitation with corresponding lapse rate corrections. Slope-dependent precipitation adjustments account for not explicitly resolved redistribution processes in mountainous terrain (Griessinger et al., 2016). Precipitation is further partitioned into rain and snow based on air temperature T_air using a logistic function (Kavetski and Kuczera, 2007; Magnusson et al., 2014; Oberrauch et al., 2024). The fraction of solid precipitation f_solid is computed as

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{f}_{\text{solid}}\left(T_{\text{air}}\right)={\left(\mathrm{1}+\exp \left({\displaystyle \frac{T_{\text{air}}-T_{\text{thresh}}}{m_{\mathrm{p}}}}\right)\right)}^{-\mathrm{1}},\end{array}$$

where T_thresh is the temperature threshold defining the transition from predominantly solid to liquid precipitation, and m_p controls the width of the mixed-phase precipitation range. The fraction of liquid precipitation is then given by $f_{\text{liquid}}=\mathrm{1}-f_{\text{solid}}$. For COSMO data, T_thresh=1.04 °C and m_p=0.15 °C are fixed parameters. Since ICON already provides solid and liquid precipitation as separate quantities, instead we infer T_thresh for each 1 km ICON grid cell such that f_solid matches the rain-snow partitioning of the ICON model output. The infered temperature threshold can then be interpolated to the 250 m model grid and applied to compute f_solid at high resolution. Note that this partitioning scheme allows for adjustments of the precipitation phase through perturbations of air temperature during the assimilation step (Sect. 2.4).

Analogous to the operational products of the OSHD, the presented dataset represents seasonal snow for individual water years, defined as the period from 1 September to 31 August of the following year. Accordingly, at the beginning of each water year, all snow is removed, disregarding the buildup of any perennial snow and firn at high elevations and on glaciers. A one-year spin-up simulation is performed to initialize the soil layer temperatures for the water year 2016. Subsequently, the end-of-season soil temperatures are used to initialize the following season.

2.4 Particle-filter-based assimilation

The particle filter (PF, Chopin and Papaspiliopoulos, 2020) is a Bayesian data assimilation method that estimates the state of a system from a weighted set of ensemble members, referred to as particles. The ensemble simulation represents the initial uncertainty of the prior estimate. Particle weights are sequentially updated based on the likelihood of the given observation, under the assumption that the respective particle represents the true state of the system. The resulting posterior distribution combines prior information and observations, accounting for their respective uncertainties. Particles are propagated from one assimilation step to the next according to the model dynamics, with a potential resampling step to ensure adequate dispersion. Given spatially correlated priors, it is possible to update unobserved locations indirectly (e.g. Cluzet et al., 2022; Odry et al., 2022; Alonso-González et al., 2023).

We employ the PF-based assimilation scheme described by Oberrauch et al. (2024) to correct spatiotemporal error patterns in the meteorological forcing data by assimilating point snow depth observations from 444 stations across the domain. Local corrections are inferred independently for each station location and three-day assimilation window, and subsequently interpolated to unobserved locations across the domain. The assimilation procedure is performed on a separate “offline” set of point simulations, allowing for computationally efficient updates of the fully distributed simulations over such a large domain without the need for a gridded ensemble (see Fig. 3 for an overview of the modeling and assimilation workflow).

Figure 3 Simplified flowchart of the modeling chain, including the PF-based assimilation scheme to correct meteorological forcing inputs.

The prior ensemble is generated by stochastically perturbing incoming longwave radiation and air temperature additively, and precipitation amount multiplicatively via a scaling factor. This perturbation strategy provides direct and largely independent handles on the radiative energy budget, as well as the amount and phase of precipitation (the latter via small perturbations of air temperature via Eq. 1), thereby minimizing equifinality issues (Oberrauch et al., 2024). Oberrauch et al. (2025) showed that the best results were achieved by allowing sufficient flexibility in correcting the given forcing variables, while additional perturbation of model parameters, namely the snow viscosity, did not notably improve the final estimate. Accordingly, to create this dataset, we applied the “METEO” perturbation strategy from Oberrauch et al. (2025), which uses less constrained perturbation priors. The additive perturbations of longwave radiation and air temperature are drawn from normal distributions with zero mean and standard deviations of ${\mathit{\sigma }}_{\text{LW}}=\mathrm{117}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and σ_TA=2.73 °C, respectively. The precipitation scaling factors are drawn from a log-normal distribution with parameters μ_PR=0 and σ_PR=1.174. In all cases, only the central 80 % of the respective distributions are sampled to exclude extreme perturbations (Oberrauch et al., 2025). Further details are given in Sect. 2.8 and Tables 2 and 3 of Oberrauch et al. (2025).

For each station and assimilation window, a set of optimal forcing corrections is defined based on the probability density distribution of the perturbation posterior. At the end of each three-day assimilation window, the posterior distribution is computed from the weight of each particle ωⁱ based on the difference between simulated snow depth $d_{\text{sim}}^i$ and observed value d_obs:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{\mathit{\omega }}^i=\exp \left(-\mathrm{0.5}{\left({\displaystyle \frac{d_{\text{sim}}^i-d_{\text{obs}}}{{\mathit{\sigma }}_{\text{obs}}}}\right)}^{\mathrm{2}}\right).\end{array}$$

The observation uncertainty σ_obs is set to 5 % of the observed snow depth, with upper and lower bounds of 5 and 20 cm, respectively. The most conservative local mode of the resulting multivariate distribution, i.e. the point with the highest probability density closest to the unperturbed state, is chosen as the optimal set of forcing corrections. For a more in-depth explanation, please refer to Oberrauch et al. (2024, Sect. 2.4.1). This deliberate collapse of the probabilistic information onto a deterministic estimate forms the basis for the subsequent propagation of the inferred information to unobserved locations.

The independently inferred local forcing corrections are interpolated spatially using a three-dimensional Gaussian interpolation scheme (Jörg-Hess et al., 2014; Oberrauch et al., 2024). The corrections at each grid point are a weighted average of all stations within a 35 km radius, weighted by the three-dimensional distance to the station locations, with vertical distances scaled by a factor γ=50 relative to the horizontal distance (following Oberrauch et al., 2024, Sect. 2.4.3). The gridded corrections obtained from the spatial interpolation are then applied to the downscaled gridded meteorological input data, forcing the distributed simulations.

The assimilation scheme notably improves snowpack simulations at unobserved locations across the domain at subregional scales, with its performance ultimately constrained by the information content of the assimilated observations (Oberrauch et al., 2025). As a result, input datasets of varying quality and processing levels from the COSMO and ICON forecasting systems can be effectively corrected, yielding snowpack simulations of comparable accuracy regardless of the input data. This homogenization of the meteorological forcing ensures consistent quality of the provided dataset throughout the entire 10-year period, preventing year-to-year discontinuities, which may arise when switching NWP models (between COSMO and ICON) or processing levels (between forecast and analysis). A detailed evaluation of the assimilation scheme, comparing multiple assimilation settings over two seasons and a range of model complexities and input data qualities, is provided in Oberrauch et al. (2025).

3 Evaluation

The OSHD model system has been successfully employed in an operational context for over a decade, delivering daily analyses and forecasts of snow cover dynamics across Switzerland, thereby providing critical information to the avalanche warning service, the Federal Office for the Environment, and other partners (Mott et al., 2023). In addition to its operational use, which relies on the system's robustness and the consistent quality of the data products delivered, continuous tuning, validation, and integration of advances from snow modeling research further contribute to the high quality of the data provided (Mazzotti et al., 2023; Quéno et al., 2024; Cluzet et al., 2024; Oberrauch et al., 2025; Haagmans et al., 2026).

To provide an estimate of the uncertainties in the presented dataset, we conducted a leave-one-station-out (LOO) validation, comparing snow depth estimates with point observations at the 444 station locations, analogous to Oberrauch et al. (2024, 2025). To ensure independence, the forcing corrections at each station location were interpolated from surrounding stations within the 35 km sphere of influence, disregarding the correction estimated at the station itself. Figure 4 shows simulated and observed snow depths across the full simulation period from September 2016 to August 2025, aggregated into four elevation bands by averaging over all stations within each band. Even in the highest elevation band, where average peak snow depths exceed 2 m, the RMSE remains as low as 6.25 cm, with a slightly negative bias of 3.34 cm.

Figure 4 Average snow depth per elevation band across all 444 station locations between September 2015–August 2025. Observed values are shown as orange dots, while model estimates from the leave-one-out validation simulation are shown as a blue line. RMSE and bias are calculated over the entire period for all snow-covered days within the respective elevation bands (i.e. excluding days when both simulated and observed values are zero).

Given the representativeness error inherent to point snow observation, a validation at individual locations is not sufficient to evaluate the dataset across the complex mountain topography. A recent validation study by Cluzet et al. (2024) compared FSM2OSHD against snow wetness observations from Sentinel-1 retrievals and identified a delayed melt onset, which was subsequently corrected through a refined albedo parameterization, now adopted in the present model configuration. Furthermore, a validation against SCF maps from Sentinel-2 retrievals presented in Oberrauch et al. (2025) showed that the PF-based assimilation scheme reliably accounts for forcing uncertainties at a subregional scale, but errors in small-scale accumulation and ablation patterns remain unresolved.

We present a spatially explicit validation of the dataset against SCF maps derived from Sentinel-2 observations over eight seasons with sufficient data availability (WY 2018 to WY 2025). For each day with a usable Sentinel-2 acquisition between March–July, we aggregate modeled and observed SCF information to produce a regional snow-line altitude (SLA) for eight aspect classes and flat locations. The SLA is defined as the elevation where the average SCF in a given region drops below 30 %. Figure 5 depicts the distribution of mean SLA bias (model-observations) over all 16 hydrological units of the Swiss Alps (as indicated in Fig. 1). Overall, the median SLA bias is within ±100 m for all water years and aspect classes, with individual subregions occasionally exceeding this range but more frequently falling well below it. Noticeable differences exist between water years, though there is no clear link to high or low snow years (cf. Fig. 6), suggesting that the bias is driven by factors other than overall snow availability. A slight aspect dependency is evident, with a tendency towards positive biases in the north-easterly sector, and negative biases in the south-westerly sector.

Figure 5 Distribution of average snow-line altitude bias over all 15 hydrological units spanning the Swiss Alps, aggregated into eight aspect classes and flat terrain. The different colors indicate the different water years (WY). Each box represents the interquartile range (IQR), with the horizontal line indicating the median value. The whiskers represent the range, while outliers, defined as points beyond 1.5×IQR from the box boundaries, appear as crosses.

Figure 6 Daily time series of average SWE (orange shading, left axis) and total meltwater runoff (blue line, right axis) over the whole domain from September 2015 to August 2025.

4 Discussion

As outlined in Sect. 3 and further detailed in Oberrauch et al. (2024, 2025), the PF assimilation scheme reliably accounts for forcing uncertainties at subregional scales, while only minor errors in small-scale accumulation and ablation patterns persist. The leave-one-out evaluation at station locations shows that snow depth estimates remain within a few centimeters of observed values across the entire 10-year period and the full elevation range. Furthermore, the evaluation against independent SCF maps demonstrates that the dataset accurately captures the heterogeneous evolution of the snowpack across the complex mountainous domain, with an average SLA bias well below 100 m in most cases. Together, the use of state-of-the-art meteorological forcing, physics-based snow modeling, and data assimilation has allowed us to produce a dataset with, to our knowledge, unprecedented spatial resolution and accuracy for a domain and period of this size. Nevertheless, a few limitations remain, which are discussed below.

Oberrauch et al. (2025) showed that assimilation performance is limited by the information content of in situ observations rather than their absolute number or density, as point observations cannot fully resolve the heterogeneous evolution of the snowpack across different slopes and aspects, especially when collected predominantly at flat-field locations (e.g. Grünewald and Lehning, 2015; Cluzet et al., 2024). Hence, to further improve model estimates, additional spatially explicit data would be necessary. Assimilating SCF observations using ensemble smoother techniques (e.g. Girotto et al., 2014; Margulis et al., 2015) is common practice in reanalysis studies (e.g. Margulis et al., 2016; Aalstad et al., 2018; Fiddes et al., 2019; Liu et al., 2021) and allows for reconstructing the seasonal snowpack evolution from observed meltout patterns. In our case, however, the available SCF data did not cover the full 10-year period, and assimilating them would have likely introduced temporal discontinuities in dataset accuracy, as reported, e.g. for the French snow reanalysis dataset (Vernay et al., 2022). Since providing a temporally homogeneous dataset has been a main objective, we chose not to assimilate SCF observations.

Snowpack estimates might be less reliable at very high elevations above 3000 m, where uncertainties are difficult to quantify due to the lack of comprehensive observational data. At these elevations, a limited number of grid cells lie outside the sphere of influence of any snow monitoring station, as defined in Sect. 2.4. For those grid cells, forcing corrections had to be interpolated from the three nearest stations, regardless of distance or elevation difference, to avoid unrealistic discontinuities in the dataset. Given that only about 1.8 % of the domain lies above 3000 m and that no alternative observational data are available for assimilation or validation at these elevations, this represents the best available approach, and any resulting errors are expected to remain limited in their impact on the overall dataset quality.

In forested areas, snow models require additional process representations and input datasets, which unavoidably introduce additional sources of uncertainty. The snow–canopy and canopy–atmosphere interactions within FSM2OSHD have been extensively validated at the process level and for hyper-resolution simulations between 1–50 m (Mazzotti et al., 2020b, 2021, 2023). At the presented 250 m resolution, a comparison with PlanetScope RGB composites has confirmed that FSM2OSHD accurately reproduces observed differences in the seasonal, interannual, and regional evolution of the snowpack among open, dense, and sparsely forested areas (Haagmans et al., 2026). To our knowledge, the extent of model validation efforts in forest terrain and the level of detail of the canopy input datasets (Webster et al., 2025) are unprecedented in the European Alps. Nevertheless, the vegetation input datasets may contain inaccuracies. For example, the vegetation data have been acquired between 2020–2022, and changes in the canopy structure since are not accounted for (Webster et al., 2025). Additionally, snow estimates in forested areas lack direct observational constraints. However, given the model's accurate representation of forest-snow processes, improvements in the meteorological forcing data will also enhance snow estimates in the forest, even when the forcing corrections were derived from stations outside the forest.

The research model variant FSM2trans (Quéno et al., 2024) includes dedicated modules for snow redistribution by gravity- and wind-driven processes (Bernhardt and Schulz, 2010; Liston et al., 2007), as well as an updated density-dependent layering scheme to represent erodible snow. While accounting for horizontal redistribution of snow enables a more realistic representation of small-scale accumulation and erosion patterns, it requires simulations at hectometre or finer spatial resolutions, along with appropriate downscaling of wind fields (Quéno et al., 2024; Berg et al., 2024; Reynolds et al., 2024). At the resolution of 250 m of the presented dataset, FSM2OSHD does not explicitly account for snow redistribution, which may affect the accuracy of local-scale snow distribution patterns.

5 Example usage

The presented dataset provides daily values of snow depth, SWE, SCF, and meltwater runoff over an area of 58 000 km² and a period of 10 years, enabling analyses across different temporal and spatial scales. Snow depth is reported as the average value per grid cell. For example, if a grid cell has a snow depth of 25 cm and an SCF of 50 %, this can be interpreted as half the cell being covered with 50 cm of snow. SWE and runoff are given as the total water volume per grid cell. Runoff represents the liquid water leaving the snowpack; during rain-on-snow events, this can include liquid precipitation percolating through the snowpack. Note that soil infiltration is not considered in these calculations (Essery, 2015; Mott et al., 2023, for details). Below, we highlight several illustrative examples of the use of this dataset.

Figure 7 Map of average pixel-wise peak SWE between September 2015–August 2025. The Swiss national border is outlined in black, and a shaded relief is used as a background map. The inset panel shows the elevation distribution of average peak SWE across 250 m elevation bands.

Figure 6 shows time series of the average SWE (orange shading, left axis) and total meltwater runoff over the entire domain (blue line, right axis). Over all 10 seasons, average peak-SWE amounts to 130 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$, with around 90 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ for low-snow years and above 160 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ for high-snow years. Interannual variations in peak-SWE values are even more pronounced, ranging from up to 197 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ in water year 2018 down to 85 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ in water year 2025. While the overall seasonality is consistent, the magnitude of seasonal accumulation and melt varies notably. On average, peak-SWE typically occurs on March 24±18 d, but can be as early as the end of February and as late as the beginning of May (corresponding to a range of 65 d). Meltwater runoff also follows a clear seasonal pattern, with low values in autumn and winter, peak values in spring, and a gradual reduction during summer. Day-to-day variability is, however, much higher, with pronounced spikes throughout the season and maximum values occurring as early as December (see water year 2021).

The average pixel-wise peak SWE over all 10 seasons, shown in Fig. 7, reveals a distinct elevation-dependent pattern. In the high mountains and glaciated areas along the Alpine Main Ridge, peak SWE values regularly exceed 1000 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ and even surpass 2000 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$. In the mid-elevation mountain ranges between 2000–3000 m, average peak SWE is approximately 550 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$. Over the Swiss Plateau, the snowpack remains shallow, with considerably higher peak-SWE values observable in the adjacent Jura Mountains.

Figure 8 Maps of snow melt-out date for the water years 2016–2025. The snow melt-out date is defined as the last day of the water year with at least 5 cm of snow depth, following a minimum of 30 consecutive days of continuous snow cover. The Swiss national border is outlined in black, and a shaded relief is used as a background map.

Figure 8 shows the snow melt-out date for each season, defined as the last day of the water year with at least 5 cm of snow depth, following a minimum of 30 consecutive days of snow cover. The Swiss Plateau rarely experiences a whole month of continuous snow cover, so a melt-out date cannot be computed for most of the region. In the Jura and Alpine foothills, the melt-out date typically falls between December–March, while snow cover can be observed up to June and later along the Alpine Main Ridge. Nevertheless, the interannual variability is considerable, and total domain-wide SWE (shown in Fig. 6) and melt-out dates are not necessarily correlated: in the low-snow year 2016, for instance, the average peak SWE was only 120 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$, yet melt-out dates above 2000 m ranged from June to September.

Figure 9 Maps of the number of snow days per season (defined as days with a minimum snow depth of 5 cm) for water years 2016–2025. Regions with 120 or more snow days and a snow cover that persists until the end of the season are labeled as >365. The Swiss national border is outlined in black, and a shaded relief is used as a background map.

Figure 9 shows maps of the number of snow days per season (defined as days with a minimum snow depth of 5 cm) grouped into discrete classes between 1, 15, 30, 60, 120, and more than 365 d. Across all 10 years, much of the Swiss Alps exhibits a persistent seasonal snowpack with more than 120 snow days, while the Alpine foothills, Jura Mountains, and Swiss Plateau generally experience far fewer snow days. The water year 2020 shows a particularly distinct bimodal pattern, with areas experiencing either more than 120 snow days or fewer than 15, and only narrow transition zones in between. The interannual variability is pronounced, not only between low and high elevations but also between the northern and southern slopes of the Alps. In most years, at least one day of snow cover occurs across the majority of the domain, although some exceptionally warm and/or dry years (e.g. 2020, 2022, and 2023) show snow-free conditions in the lowlands. Regions with 120 or more snow days and a snow cover that persists until the end of the season (labeled as >365 in the figure) are confined to small areas in the highest and glaciated parts of the Alps and are virtually absent in certain years, such as 2017.

6 Data format

The dataset provides daily estimates of snow depth (m), snow water equivalent ($\mathrm{mm}\mathrm{w}.\mathrm{e}.\equiv {\mathrm{kgm}}^{-\mathrm{2}}$), snow cover fraction (unitless), and snowmelt runoff ($\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}\equiv \mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) for each grid cell within the domain. The reference time for daily values is 06:00 Central European Time ($\text{CET}\equiv \text{UTC}+\mathrm{0}$1:00), with snow depth, SWE, and SCF representing states at that time, and runoff corresponding to the accumulated total over the preceding 24 h. Given that on 1 September at 06:00 CET of each season, all snow is removed, all variables for that day in the dataset are zero because the end-states from the previous season are not carried over.

The data are stored in self-explanatory monthly NetCDF files (e.g. OSHD_DATA_2020-01.nc for January 2020) within individual zip archives for each water year (e.g. OSHD_DATA_WY2020.zip for water year 2020), with a total data volume of about 800 MB per season. Variable names, attributes, and metadata adhere to the CF 1.12 and ACDD 1.3 conventions. The data variables are structured as three-dimensional arrays, with time, easting, and northing coordinates. The temporal coordinate denotes the number of days since the first day of the corresponding month. The horizontal coordinates refer to the center of the respective grid cell, given in the local Swiss CH1903+/LV95 reference system (EPSG:2056). The total spatial extent is also specified by the minimum and maximum longitude and latitude in WGS84 coordinates (EPSG:4326). Additional metadata includes the dataset title, keywords, summary, and version history; contact information; and a reference to the FSM model version used (as GitHub tag). The model DEM is provided as a GeoTiff (OSHD_MODEL_DEM.tif). Compatibility with the NASA Panoply NetCDF Viewer (NASA GISS, 2025) was verified using Panoply v5.7.1 under Windows 11.

7 Conclusions

We present a reanalysis dataset of spatially explicit seasonal snow cover dynamics for Switzerland and its bordering regions for water years 2016–2025, based on the high-resolution simulations with the physics-based, multi-layer snow model FSM2OSHD and the assimilation of snow depth observations from 444 stations across the domain. The combination of dedicated downscaling of NWP forcing data with the upscaling of state-of-the-art hyper-resolution DEM and light-availability datasets enables simulations that account for subgrid variability within the 250 m model resolution, which is crucial for capturing the spatial heterogeneity of the mountainous snowpack. The coupled energy- and mass-balance equation, solved for individual numerical snow layers, explicitly accounts for different atmospheric and snowpack processes of open, forested, and glaciated areas. Finally, the assimilation of in situ snow depth observations dynamically corrects spatiotemporal error patterns in the meteorological input data, ensuring consistent forcing quality over the presented 10-year period.

Spanning 10 years of physically consistent snow cover estimates across a large part of the European Alps, the dataset contributes to the snow data record for a region characterized by complex terrain and diverse hydroclimatic regimes. With its high spatial and temporal resolution, the dataset helps bridge the gap between coarse global products and detailed local observations. The data are publicly and freely available at https://doi.org/10.5281/zenodo.17313889 (Oberrauch et al., 2026), providing a robust foundation for hydrological, ecological, and cryospheric analyses.