The physics-based FSM2oshd model (https://github.com/oshd-slf/FSM2oshd) solves the complete mass and energy balance of the snowpack for open and forested areas (Mazzotti et al., 2020a; Mazzotti et al., 2020b) at an hourly temporal resolution (Figure 1, FSM2oshd). FSM2oshd was originally based on the FSM2 snow cover model (Flexible snow model, https://github.com/RichardEssery/FSM2), with adapted snowpack process parameterizations and new model components for snow cover fraction (Helbig et al., 2021a) and for forest processes (Mazzotti et al., 2021; Mazzotti et al., 2023; Webster et al., 2023). In its operational setting the physics based model chain is run at a horizontal resolution of 250 m, but is also applied at higher resolutions (up to 25 m) for research purposes. The current physics-based OSHD model chain was initiated by Magnusson et al. (2015) who built a Swiss snow model framework suitable for predicting daily SWE and snowmelt using the first multi-model framework of physically based energy-balance models, JULES Investigation model (JIM), from Essery et al. (2013). The spatialized version of JIM was set up to run with downscaled numerical weather prediction data (NWP) and became the first operational OSHD version.

The performance of JIMoshd was evaluated in a study by Griessinger et al. (2019), who ran the latest version of the model at a spatial resolution of 250 m and temporal resolution of 1 h. The input was then fed into hydrological StreamFlow model simulations for 25 catchments in Switzerland. The study found that daily discharge simulations using StreamFlow were more accurate when JIMoshd simulations were updated with snow model mass and energy fluxes using snow observations with optimal interpolation (OI; Magnusson et al., 2014) or the bias-detecting ensemble approach (BDE; Winstral et al., 2019). In another evaluation, JIMoshd simulated snow cover evolution at a spatial resolution of 1 km and showed good seasonal agreement with various spatial snow cover data. This ranged from airborne-acquired fine-scale snow depth data to satellite and terrestrial imagery for Switzerland (Helbig et al., 2021b).

Essery (2015) developed the initial version of the Factorial Snow Model (FSM) to provide a more organized approach to the development of the snow model suite JIM. To enhance the model’s physics, particularly adding the representation of forest snow processes (for more details see Section 2.1.3) and output options, Essery (2015) introduced the Flexible Snow Model FSM2. The FSMoshd, originally based on FSM2, is capable of running on both a list of points and a regular grid simultaneously. Grid cells are divided into separate tiles, including forest-covered, open, and glacierized fractions (Figure 2, left), with fractions covered by large bodies of water where snow does not accumulate being excluded. The states and fluxes of each tile are computed as weighted averages of the respective states and fluxes of the grid cell. FSM2oshd includes local canopy cover fraction, hemispherical sky view fraction, wind attenuation, spatial patterns of shortwave and longwave radiation transfer through the canopy, and preferential deposition of snow in canopy gaps for the forest tiles, which have been improved by Mazzotti et al. (2020a), Mazzotti et al. (2020b). Station locations are linked to simulations at the point scale to assimilate local observations into the model (see Section 3.1).

The FSM2oshd model, which is based on physics, is currently being used for seasonal runs of recent water years (since WY 2016) to analyze and predict the spatio-temporal dynamics of snow water equivalent (SWE), snow depth (HS), snow melt runoff (SMR), and snow state (based on liquid water content) on a daily and hourly basis (Figure 3). Additionally, a new snow height model, FSM2oshd-HN, has been implemented to explicitly model the height of new snow over defined time frames, taking into account snowpack processes such as melting and compaction (Figure 4). The HN model simulates a HN board, where snow accumulates and settles over a specific time period. At the beginning of the accumulation period, the bottom boundary condition of the snowpack is initialized based on snow surface temperatures at the initialization times. In forecast mode, the HN model runs for lead times of 12, 24, 48, and 72 h, driven by COSMO forecasts (Figure 4). Due to the settling of the snowpack, the HN may be smaller after 48 h than after 24 h, despite continued snowfall. Furthermore, the HN model is run in the hindcast, using the data assimilation scheme described in Sect three. To allow for an optimal transition between hindcast and forecast data, precipitation is obtained from the nowcasting product INCA (Integrated Nowcasting through Comprehensive Analysis) (Meteoswiss, 2020) to drive the FSM2oshd-HN model. All other meteorological variables are forced by the COSMO analysis or forecast models. Within the operational model framework, the HN model is updated whenever new forecast and analysis data from COSMO becomes available.

The physics-based OSHD model chain has a wide range of applications, including operational use as well as in various research projects. These projects aim to study high-resolution snow melt processes in alpine terrain (Griessinger et al., 2019) and forested areas (Mazzotti et al., 2023), and to project future snow water resources and rain-on-snow events (e.g., Arnoux et al., 2021; Schirmer et al., 2022). The advantage of using a physics-based model is that it allows for more precise predictions of snow water inputs during snow melt events, particularly those driven by turbulent heat fluxes such as rain-on-snow situations. Additionally, the model can resolve radiation differences between slopes, which significantly increases the spatial variability of snow melt in the spring season (see Figure 5). Compared to temperature-index based models, which typically predict much lower snow melt during rain-on-snow events (as shown in Figure 6), the physics-based model chain appears to better capture snow melt processes during complex meteorological conditions.

The physics-based model chain used for snowpack simulation relies heavily on precise and high-resolution information about the local meteorological conditions that affect the mass and energy balance of the snowpack, both spatially and temporally. For instance, the FSM2oshd model requires hourly data on various parameters such as air temperature, air pressure, relative humidity, total precipitation, wind speed, shortwave and longwave radiation, for each grid cell (as shown in Figure 1). Although such detailed information can be obtained by downscaling the high-resolution (1–2 km) COSMO forcings, it is important to note that these forcings may still contain errors inherited from the atmospheric models. This is especially true for mountainous regions, where NWP models often lack assimilation of surface observations.

To minimize errors in the simulation, snow data assimilation can be used to improve the accuracy of solid precipitation by correcting it before running the snowpack model (front-end) or during the execution of FSMoshd (back-end), utilizing the wealth of snow depth and 2 m air temperature data available in Switzerland and neighboring countries. The procedure for solid precipitation data assimilation is detailed in Section 3 (as shown in Figure 1).

The conceptual model chain was developed within the OSHD model chain (Figure 1, OSHD_TICL) for climatological runs to compare snow water resources in Switzerland of water years since 1998 at a horizontal resolution of 1,000 m (Figure 5). The conceptual model is based on an efficient temperature-index (TI) snow model with SWE and liquid water as main state variables and only requires precipitation and temperature as forcing data, that is, available for a longer historical period than the physics-based model chain. The model is described in detail in Magnusson et al. (2014); Griessinger et al. (2016). The conceptual models include a density model (HS-2-SWE, see Section 3.1.1) for converting observed changes in snow depth to solid precipitation fluxes and considering decreasing snow depth due to compaction (Jonas et al., 2009). The temperature index method is used to compute the daily amount of snow melt occurring.

Our model framework uses an Ensemble Kalman filter (EnKF) to update simulated melt rates and two model state variables (SWE and liquid water content) based on observed snow water equivalents computed from snow depth records. We run 100 ensemble members for the point observation sites and the model grid, using interpolated input data perturbed with spatially correlated errors (Magnusson et al., 2014).

The model framework employs two conceptual snow models, OSHD_TICL (Climatological Temperature Index Model) and OSHD_ETI (Enhanced Temperature Index Model). In our operational setting, we use OSHD_TICL, which compares SWE of different years using the same data sources available since 1998 (Figure 7). We ensure maximum temporal constancy between individual water years since 1998 by considering only stations and data sources that feature data since 1998. Thus, OSHD_TICL only uses gridded daily precipitation data (RHIRESD; MeteoSwiss, 2019) and observed solid precipitation fluxes. Further details on data assimilation procedures for precipitation are provided in Section 3. To account for the sub-grid spatial variability of snow, we use the seasonal and scale-independent algorithm for the fraction of snow-covered area (fSCA) developed by Helbig et al. (2015) (see Section 2.1.2).

We also developed the non-operational ETI model to improve the representation of the precipitation phase. This model uses hourly air temperature and precipitation data to discriminate between rain and snowfall. Forcing data for air temperature, relative humidity, and shortwave radiation is downscaled from respective fields from COSMO (see Section 4).

A full data set for climatological runs, providing maps of snow water equivalent and snow cover fraction from water year 1998–2022 from the OSHD_TICL model, is available at EnviDat. doi:10.16904/envidat.401 .

In order to gain information on daily solid precipitation rates, we convert daily snow depth records (HS) to snow water equivalent (SWE) using a density model (HS-2-SWE) (step 1, Figure 1). This method uses station data (Figure 8) from IMIS (Intercantonal Measurement and Information System) snow monitoring networks including automatic stations and manual observer data or from SwissMetNet (MeteoSwiss). Furthermore, data from DWD (Deutscher Wetterdienst, German national meteorological service), ZAMG (Zentralanstalt für Meteorologie und Geodynamik, Austrian national meteorological service) and MeteoFrance (French national meteorological service) are operationally used for regions outside the Swiss border. The station network consists of 444 stations measuring snow depths (HS stations) and 256 stations measuring air temperature (TA-stations). For the climatological runs the number of HS stations is reduced to 349 stations and TA-stations to 255 providing a full time-series of quality-controlled data from WY 1998 until today. Snow depth and air temperature data are extracted at daily timesteps and are quality checked.

We convert daily point measurements of snow depth to water equivalent using a model based on the methods presented by Martinec and Rango (1991). The adapted model describes the accumulation and settling of the snowpack for each layer. The density model applied for each layer has been calibrated based on data from the over 10,000 snow profiles used in the study by Jonas et al. (2009). We can compute snowfall and melt rates from the increments/decrements in snow water equivalent. For an example of the computed solid precipitation and melt time series, see Figure 2 in Magnusson et al. (2014). Note that melting and accumulation cannot occur simultaneously with this simplified description of the snowpack development. At the same time, this model only requires daily snow depth records without the need for any meteorological data.

To carry out the assimilation outlined in Step 3 of Figure 1, obtaining the best possible prior estimate for solid precipitation on the model grid is crucial (Step 2, Figure 1). To achieve this, we gather data from various sources, including gridded daily precipitation products (RhiresD) from MeteoSwiss and interpolated temperature fields based on large station datasets. RhiresD provides a spatial analysis of daily precipitation over the OSHD domain, utilizing data from heated rain gauges, typically 430 within and 220 outside of Switzerland daily. The daily temperature interpolation follows the methods presented by Frei (2014) and employs data from over 300 stations covering most of the elevation range given by the OSHD domain. We combine these daily products with hourly data from COSMO 1 E to improve, for instance, precipitation phase determination to obtain the best possible priors of daily snowfall on the OSHD grid, which are further used in the optimal interpolation step described below.

The gauge network utilized in the gridded precipitation products of MeteoSwiss has a limited representation of elevations above 1,200 m. Additional snow depth data from the snow observation networks, particularly the IMIS network, are incorporated to enhance the accuracy of solid precipitation information at high altitudes. The first-guess or background fields of daily snowfall obtained in the previous step (refer to Section 3.2) serve as the basis for assimilating point solid precipitation rates inferred from snow depth measurements (refer to Section 4.1). An optimal interpolation scheme (OI) is employed for the assimilation, efficiently updating the background field utilizing observations (step3, Figure 1). The OI method considers apriori-defined covariances between errors in the gridded input fields and point observations, accounting for observation errors. We assume errors in the precipitation input fields to be spatially correlated, while observation errors are assumed to be spatially uncorrelated (Magnusson et al., 2014). The optimal interpolation scheme generates a best-guess daily solid precipitation field, incorporating information from the snow depth stations.

To account for the effects of terrain on precipitation and snow redistribution by wind and avalanches, distributed precipitation multipliers are utilized as snow observations are commonly conducted at flat field sites. To implement the sub-grid precipitation adjustment, snow depths obtained from airborne lidar acquisitions in the European Alps, as presented in Grünewald and Lehning (2015), are employed. The scaling method, following the iterative scaling approach proposed by Vögeli et al. (2016), corrects for snow accumulation and is based on slope calculations conducted at 25 m and then upscaled to a resolution of 250 m. It should be noted that the scaling method is only applied to solid precipitation.

In Figure 9, the impact of the assimilation procedure described above on the accuracy of snow simulations at 45 measurement locations spanning an elevation range of 1,195–2,690 m.a.s.l. from 2015-09-02 to 2022-8-01 is presented. The optimal interpolation step leads to a reduction in the root-mean-squared error and bias for all elevation bands. It should be noted that the error reduction is smallest for the two highest altitude bands, likely due to lower station density at these elevations and, therefore, fewer data available for assimilation. In Figure 10, we present a comparison between the time series of modelled snow depths with and without applying OI technique, based on observations for Water Year (WY) 2022. To ensure robust analysis, station data is aggregated for 500 m elevation bands. The results clearly highlight the significant benefits of employing optimal interpolation in operational snow modeling, particularly in reducing errors associated with solid precipitation. These findings further corroborate and validate the earlier conclusions reported in studies conducted by Magnusson et al. (2014) and Griessinger et al. (2019).

The OSHD-EB model chain employs the Bias-Detecting Ensemble (BDE) method to enhance snow accumulation and depletion predictions by updating mass and energy fluxes with snow observations (Winstral et al., 2019). The BDE method evaluates the entire modeling chain to identify potential biases in the meteorological driving data and snow model, serving as an indirect approach for detecting biases in solid precipitation and surface energy inputs to the model.

In the first step, the BDE algorithm determines whether accumulation, melting or settling dominates the observed snow depth development. An ensemble of simulations is then conducted and evaluated against the snow observations, with precipitation perturbations evaluated on snow days and only energy flux perturbations on melt days. Daily history of the best-performing biases based on simulated and measured snow depths at each site is updated, with the assumption that these biases exhibit temporal coherence. Daily values are back-averaged and spatially interpolated to the modeling grid (Step 4, Figure 1).

Griessinger et al. (2019) have demonstrated that the use of the BDE method significantly improves the model performance of runoff modeling results for energy balance snow simulations.

For model applications with target resolutions at the sub-kilometer scale, radiation is dynamically downscaled using topographical information from a 25-m resolution digital elevation model accounting for topographic shading and the inclination of the respective grid cell following Jonas et al. (2020). To this end, the DEM is evaluated for each location of interest and point in time to determine whether the Sun is visible (weather permitting) or obscured by terrain. This enables calculating incoming shortwave radiation within terrain from either the diffuse component or the diffuse and direct components of the above-terrain radiation. The local hemispherical sky view fraction is used to constrain diffuse shortwave radiation and partition incoming longwave radiation into the fractions sourced from the sky and the terrain. In the operational setting, radiation per inclined surface is re-upscaled from 25 m to the current target model resolution of 250 m (Figure 11).

For higher-resolution model applications (sub-kilometre scale), two wind downscaling methods are available within the OSHD model framework. The two statistical approaches (Helbig et al., 2017; Winstral et al., 2017) aim to downscale gridded coarse-scale wind data to finer grids using information from the local topography. However, these two methods differ in the required wind data and topographic parameters. In the current operational model framework, the approach suggested by Winstral et al. (2017) is used.

The statistical approach, Winstral et al. (2017) developed uses local topographic information at a much higher resolution than present in the forecast wind data to delineate the spatial organization of biases in the forecast products. The local terrain is statistically linked with observed forecast biases at stations located at different terrain shading locations. The sub-grid or high-resolution terrain information is determined from a 25 m resolution DEM using two established terrain parameters: 1) the topographic position indices (TPI), which describes elevation relative to its surroundings, independent from wind direction (Jenness 2006); 2) the terrain parameter Sx which is a direction-dependent, slope-based assessment of topographic shelter and exposure capable of differentiating such slopes based on given wind directions. The Sx parameter has been used successfully to distribute hourly observed wind speeds across a modeling domain (Winstral and Marks, 2002; Winstral et al., 2009; Winstral et al., 2013; Fiddes and Gruber, 2014). The resulting wind exposure metrics effectively segregated high (ridges), moderate (slopes), and low wind speed sites (valleys). An optimization scheme is applied to further improve downscaled wind distributions, nudging the coarse-scale wind distributions from COSMO data to better match the observed wind speeds at stations (distributions and variability). The optimization method requires a reasonable number of wind stations in the domain and a new optimization whenever new forecast data is used.

As an alternative approach, the statistical downscaling scheme, described by Helbig et al. (2017), is available in the OSHD model framework. This method can downscale coarse-scale wind speeds from numerical weather prediction models to points or high-resolution grids by describing the impact of the surrounding resolved topography. Easily derived terrain parameters using the mean squared slope for each fine-scale grid and a parameter related to the Laplacian of terrain elevations are used to calculate a local topographic downscaling factor. Both terrain parameters relate to the sheltering and exposure of grid cells. Depending on the two terrain parameters, the downscaling scheme is based on millions of simulated ARPS (Advanced Regional Prediction System) fine-scale wind speed values on model topographies (Gaussian random fields), depicting a wide range of realistic terrain characteristics.