Application of physical snowpack models in support of operational avalanche hazard forecasting: A status report on current implementations and prospects for the future

The application of numerical modelling of the snowpack in support of avalanche hazard prediction is increasing. Modelling, in complement to direct observations and weather forecasting, provides information otherwise un-available on the present and future state of the snowpack and its mechanical stability. However, there is often a perceived mismatch between the capabilities of modelling tools developed by research organizations and implemented by some operational services, and the actual operational use of those by avalanche forecasters. This causes frustration on both sides. By summarizing currently implemented modelling tools specifically designed for avalanche forecasting, we intend to diminish and contribute to bridging this gap. We highlight specific features and potential added value, as well as challenges preventing a more widespread use of these modelling tools. Lessons learned from currently used methods are explored and provided, as well as prospects for the future, including a list of the most critical issues to be addressed.


Introduction
Avalanche hazard forecasting requires information about the past, current and future state of the snowpack. This includes accounting for the vertical profile of its key microstructural and mechanical properties. Avalanche hazard forecasters have traditionally focused on field observations of snow conditions, despite the significant spatial variability of snow conditions and the potential non-representativeness of point observations. Such observations are based on observations from registered observers at established observation station and field reports, automated observation networks, outing reports on mountaineering and skitouring community websites and observations from the forecasters themselves. The extrapolation in space and time is typically based on the forecaster's ability to combine the broad diversity of snow and meteorological observations and forecasts, together with knowledge relevant to snow evolution processes. Physically-based snowpack models have been developed since the 1980s, in order to to provide avalanche forecasters with information complementary to field observations and meteorological forecasts. Several physically-based snowpack models initially dedicated to avalanche forecasting purposes https://doi.org/10.1016/j.coldregions.2019.102910 Received 30 March 2019; Received in revised form 5 July 2019; Accepted 1 October 2019 have been developed for the last decades, such as Crocus originally in France (Brun et al., 1989(Brun et al., , 1992Vionnet et al., 2012) and SNOWPACK originally in Switzerland (Lehning et al., 1999;Bartelt and Lehning, 2002;Lehning et al., 2002aLehning et al., , 2002b. The physical principles upon which they are based are rather close. They however differ in the way they have been implemented for operational activities in their host and collaborating organizations. Differences are found in terms of nature and use of their meteorological driving data (e.g. balance between point-scale meteorological observations and output of numerical weather prediction -NWP-models), but also regarding the way model outputs have been post-processed and provided to the operational workstation of avalanche forecasters.
It has become apparent in recent years that an increased number of avalanche forecasting services are considering using physically-based snowpack models in support of their operational activities (e.g. Vikhamar-Schuler et al., 2011;Floyer et al., 2016). However, there has been no assessment hitherto on the successes and lessons learned regarding the use of such models for operational applications. Indeed, scientific publications tend to focus on the description of newly developed model chains, and mostly address quantitative assessments of their predictive performances against meteorological or snow observations, and not necessarily regarding their added-value for avalanche forecasting itself.
The present article contributes to bridging the gap between the research community, which has devoted significant efforts to the development of snowpack modelling chains and proposed visualization of output data, and operational avalanche forecasting centers, which have gathered experience and expressed challenges about the use of such models. This article describes the snowpack models used operationally, their operational setup, introduces some examples of their operational use and shed some light on challenges associated with using these models in an operational forecast environment, using an Information Quality framework (Bovee et al., 2003). In addition to the detailed snowpack models Crocus and SNOWPACK, we also describe simpler modelling tools SNOWGRID and seNorge, which are also used by avalanche forecasting centers in Austria and Norway, respectively.

Existing one-dimensional snowpack models used operationally
Snow on the ground evolves constantly due to exchange processes at its boundary with the overlying atmosphere and underlying ground, and under the action of internal transformation processes referred to as snow metamorphism (Armstrong and Brun, 2008). Interfacial energy and mass fluxes and internal processes are strongly coupled. For example, the surface energy balance is significantly influenced by snow albedo, which depends on near surface snow impurity content and microstructure properties. The latter depend on the physical evolution of the snow and in particular the temperature field near the surface, which in turn depends on the surface energy and mass balance. Seasonal snow on the ground can remain there for several months, so that its status at one point in time may depend on the chronology of meteorological conditions and their interaction with snow processes several weeks to months before. Therefore, meteorologically driven snowpack modelling needs to either predict the unfolding of snow conditions starting from snow-free ground at the beginning, or the capability to initialize snow conditions at a given time, which requires the capacity to measure all of the model prognostic variables at the start date of the numerical simulation. Adequately representing in a numerical model the main energy and mass fluxes at the snowpack interfaces, assuming planar layering geometry, requires several physical ingredients such as those able to capture the variations of snow density and albedo, and account for the internal energy storage associated to phase change processes in the snow (Essery et al., 2013). Given the significant vertical variations of physical snow properties in most observed snowpacks, and the space-time coupling of snow processes such as heat conduction driven by diurnally variable atmospheric boundary conditions, a multi-layer approach is generally considered appropriate for representing snow in a physically-based numerical model (Colbeck, 1991;Armstrong and Brun, 2008).
Furthermore, the likelihood of avalanche release depends on the mechanical stability of the snowpack which is closely linked to the layering characteristics of snow properties (snow stratigraphy). In this context, structural weaknesses within the snowpack are of particular importance. These are typically characterized by specific microstructural contrasts within the snowpack, such as large differences in grain size and hardness (van Herwijnen and Jamieson, 2007;. This is particularly true for structures involving a cohesive snow slab overlying a so-called weak layer. The latter is typically several mm to cm thick and made of fragile and coarse grained snow (typically surface hoar, facetted crystals or depth hoar) and originate from various processes (surface energy balance for surface hoar, temperature gradient metamorphism for facetted crystals and depth hoar) . Thus, a physically-based snowpack model for avalanche hazard forecasting must be able to accurately reproduce the most significant processes responsible for the inception of mechanically unstable conditions, which requires the capability to predict the development of thin weak layers within the snowpack.
What has become referred to as a detailed physically-based snowpack model used in support of avalanche hazard forecasting is thus a numerical model using a potentially large (i.e., more than 20) number of numerical layers for the representation of the vertical profile of snow properties, including prognostic representations of thermal (temperature, liquid water content) and microstructural (density, snow type descriptors) state of each of the layers, and equations representing their (coupled) time evolution. Simpler models, using a smaller number of layers and therefore representing a smaller number of internal processes, are also used in avalanche forecasting centers, although the breadth of model outputs they provide is more limited, and concerns mostly snow height and mass changes, and surface temperature evolution.
Regardless of the level of complexity of the 1D snowpack model employed, they cannot capture some key processes shaping the mountain snow cover. Lateral snow redistribution processes, driven by wind erosion and re-deposition, exert a strong influence on snow conditions (Mott et al., 2010;Vionnet et al., 2014;Helbig et al., 2017). They induce deviations from planar layering geometry. Local topography and interactions with the vegetation also induce such deviations. Physically-based models have been developed to represent explicitly the interaction with topography or vegetation (Pomeroy et al., 1993;Durand et al., 2005;Mott et al., 2010;Vionnet et al., 2014;Krinner et al., 2018), however, none are known to be currently used operationally by avalanche forecasting centers. It has been considered, hitherto, that there is a lower need to account for snow/vegetation interactions for avalanche warning purposes, largely due to the fact that avalanche release areas are generally located in non-forested areas. Only one-dimensional models assuming layering parallel to the local slope and operating outside forested areas, such as Crocus, SNOWPACK, SNOWGRID and seNorge are used operationally, and are described below with more detail.

Crocus
The detailed snowpack model Crocus was developed at Météo-France in the late 1980s (Brun et al., 1989(Brun et al., , 1992. Although this was never used operationally in this context, the initial purpose of this model was to interpolate in time weekly snowpit observations. Building on previously developed numerical models for seasonal snow on the ground (Anderson, 1968(Anderson, , 1976Navarre, 1975), the main modelling principles behind Crocus are given below: • The model solves the heat diffusion equation in a stratified snowpack. The thermal conductivity of snow is parameterized as a S. Morin, et al. Cold Regions Science and Technology 170 (2020) 2 function of density following Yen (1981).
• Layers are handled in a Lagrangian way, i.e. one or several snow layer(s) is (are) created upon a snow precipitation event. However, to maintain a reasonable number of numerical layers, layers can be merged if they feature sufficiently close physical properties. This dynamical layering mechanism also maintains a higher vertical resolution near the surface and at the bottom of the snowpack so as to place emphasis on locations with large vertical gradients, and accurate representation of critical physical processes occurring there (Brun et al., 1989(Brun et al., , 1992Vionnet et al., 2012).
• Snow type is diagnosed from variables representing snow microstructure (specific surface area, sphericity, and a "history" variable), whose time evolution is parameterized as a function of temperature, temperature gradient and liquid water content (Brun et al., 1992;Vionnet et al., 2012;Carmagnola et al., 2014).
• The model is purely driven by meteorological conditions at a sufficient time resolution (typically hourly). This concerns air temperature, specific humidity, wind speed, snow precipitation amount, liquid precipitation amount, incoming shortwave radiation, incoming longwave radiation.
• Fresh snow physical properties are parameterized as a function of air temperature and wind speed (Vionnet et al., 2012;Lafaysse et al., 2017).
• Snow albedo is diagnosed by the model for three visible/near infrared spectral bands depending on near-surface snow properties (age of the layer as a proxy to impurity content, optical radius computed from the snow microstructure variables) (Vionnet et al., 2012).
• All surface energy and mass balance terms are computed by the model, including surface snow sublimation or atmospheric water vapor condensation associated with latent heat exhanges. However, in the latter case which is generally associated with surface hoar formation, mass is simply added to the uppermost snow layer without changing its physical properties (Vionnet et al., 2012).
• Water percolation is based on a bucket approach defining a maximum liquid water holding capacity (Vionnet et al., 2012;Lafaysse et al., 2017).
• Compaction is based on a visco-elastic model with an empirical formulation of the viscosity (Vionnet et al., 2012;Lafaysse et al., 2017).
Since 2015, the original Crocus model has been replaced for operational use by a newer implementation of the model (Vionnet et al., 2012;Carmagnola et al., 2013;Lafaysse et al., 2017), which is now fully embedded in the land surface model ISBA within the SURFEX interface (Masson et al., 2013). One of the main differences between these two model implementations is the fact that Crocus is now fully coupled to a physically-based soil diffusion scheme (Decharme et al., 2011). Before that, the boundary condition at the bottom of the snowpack was based on a climatological heat flux, which did not depend on the current snow and meteorological conditions. The newer version allows a sub-freezing soil/snow interface and is expected to better account for the impact of snow-free conditions before the inception of the seasonal snowpack. Durand et al. (2001) developed a representation of snow drifting (erosion/submimation/redeposition) applicable to idealized geometries. In their model, referred to as SYTRON, slopes with opposing aspect but similar elevation and slope angle can exchange snow mass depending on wind speed and direction, and on simulated properties of snow on the ground from Crocus (driftability index, see Vionnet et al. (2013) for details). It has been running in real time and provided to the forecasters since the winter 2015-2016 .
Crocus (Brun et al., 1989(Brun et al., , 1992 and SURFEX/ISBA-Crocus (Vionnet et al., 2012) are both written in FORTRAN90, consistent with the widespread use of FORTRAN in most NWP models, which inspired several of their components and with which they can be coupled. MPI instructions allow for parallelization for large domains applications.

SNOWPACK
SNOWPACK has been designed in the late 1990s at SLF Davos by Lehning et al. (1999), following previous snow modelling developments such as Daisy (Bader and Weilenmann, 1992) and the use of Crocus (Fierz et al., 1997) for Swiss locations. Initially written in C, SNOWP-ACK is now a C++ stand-alone, open source application  that can also be used as a library, for example, with the Alpine3D land surface scheme .
An introduction to the basic principle of SNOWPACK can be found in the model forge models.slf.ch and a good overview of the model for research applications is given in Wever et al. (2014). Below we summarize the key features of the model with regard to its operational implementation: • SNOWPACK is a multi-layer, one-dimensional physically-based model, which solves the instationary heat transfer equations using a Lagrangian finite element method. It does not impose a fixed limit in the number of numerical layers and can include soil layers as well (see for example (Wever et al., 2014)).
• Snow is modelled as a three-phase ice, water and air porous medium, characterized by the volumetric contents of each phase. SNOWPACK shares some basic principles with the original Crocus version (e.g. the representation of snow microstructure using dendricity, sphericity and grain size, and the general form of snow metamorphism evolution laws).
• New snow density is parameterized as a function of air temperature, relative air humidity, and wind speed (Schmucki et al., 2014).
• Shortwave radiation is modelled as a volumetric heat source. The intensity of the absorbed shortwave radiation decreases exponentially with distance below the surface. Albedo is parameterized as described in Lehning et al. (2002b).
• The explicit representation of surface hoar growth (Lehning et al., 2002b), leading to the creation of additional snow layers if surface hoar formation exceeds a given threshold before a snowfall.
• The explicit erosion of snow by wind if measurements of wind speed and snow height are available (Lehning and Fierz, 2008) In general SNOWPACK uses input at an half-hourly time step and the MeteoIO library prepares the driving data sets (Bavay and Egger, 2014), including filtering of erroneous data and, for missing data, temporal interpolation or gap-filling by spatially interpolating from neighbour data sources. To drive the accumulation of snow, SNOWP-ACK can be run using either a time series of observed snow height or a time series of precipitation. In general, solid precipitation is assumed if air temperature is below a given threshold, typically 1.2°C. Regarding surface boundary conditions, SNOWPACK can use either a mix of Dirichlet and Neumann conditions or Neumann conditions only. In the former case, SNOWPACK switches from prescribed snow surface temperature (Dirichlet) to prescribed energy fluxes at the surface (Neumann) whenever the snow surface temperature exceeds a threshold, typically −1°C. Net shortwave radiation, however, is required in both configurations. Either parameterized or measured albedo can be used to estimate net short wave radiation.
The required parameters are summarized in Table 1 for different set-ups used for operational use. Here we provide more information on the virtual slope column (last column on the right). Precipitation amount as well as direct and diffuse incoming short-wave radiation similar to flat field simulations are used, and precipitation and direct short-wave radiation are projected on the slope. In addition, incoming long-wave radiation can be estimated from the closure of the energy balance on the flat field. Together with air temperature, relative humidity and wind speed, simulations on so-called virtual slopes are performed for each station. Operationally, SLF uses four aspects (N, E, S, W) at a slope angle of 38°at the same elevation as the measurement station. These virtual slopes are then exposed to wind measured on nearby ridges to estimate the amount of snow transported by the wind during a 24-hour period (Lehning and Fierz, 2008). Virtual slope simulations provide the same type of output as the simulations in the flat terrain.

SNOWGRID
The physically-based and spatially distributed snow cover model SNOWGRID (Olefs et al., 2013) is mainly developed at Zentralanstalt für Meteorologie und Geodynamik (ZAMG, Austria) and is based on a simple 2-layer scheme, considering settling, the heat and liquid water content of the snow cover, snowline depression effects and the energy added by rain. For every time step and layer, the state variables snow density, snow water equivalent, snow temperature, liquid water content, bottom liquid water flux and the surface albedo are calculated and stored. The primary focus of the model is to provide fast calculations on a large grid (see section 2.3.3) and to accurately represent the spatial distribution of these snow variables. At present, a detailed multi-layer description of the snow cover is not foreseen to be implemented in SNOWGRID, and neither are snow stability indicators based on the microstructure of the snow cover.
SNOWGRID can be operated in three different complexities: a simple degree-day scheme (energy balance parameterized by the air temperature), a more complex extended degree-day model (adding global radiation on the real surface and surface albedo; following Pellicciotti et al., 2005) and a full energy balance mode considering all components of the energy balance (the latter is still in development). The two more complex modes comprise partly newly developed schemes (e.g. for radiation and cloudiness) based on high quality solar and terrestrial radiation data (Olefs et al., 2013, satellite products and ground measurements.

seNorge
The seNorge snow model component is a single layer snowpack model which simulates different snow-related variables using a degreeday approach, such as snow water equivalent, 24 h and 72 h new snow amounts, snow height, bulk snow density and the amount of liquid water in snowpack (Saloranta, 2012). It requires either daily or 3 h (high resolution mode) mean air temperature and sum of precipitation. The model consists of two main modules, a SWE (Snow Water Equivalent) model for snowpack water balance and a snow compaction and density module used to convert SWE to SD (Snow Depth). The SWE model for snowpack water balance is run before the snow compaction and density module (Saloranta, 2012). Each output is available on a 1 km grid covering Norway and large parts of Sweden at a temporal resolution of either 3 h or 24 h. Maps of new snow amounts are an important part in the daily assessment of avalanche danger for the Norwegian Avalanche Warning Service.

Currently operationally implemented meteorological forcing configurations and associated geometry
Physically-based snowpack models operate intrinsically at the point scale, i.e. they are driven by a number of individual time series of colocated meteorological variables at the sub-diurnal time resolution. They predict the time evolution of the physical properties of snow corresponding to these time series. Several approaches have been developed to generate such meteorological driving data. Snowpack modelling for avalanche hazard prediction can be considered a demanding subclass of hydrological modelling in mountain regions. However, snowpack modelling for avalanche hazard application generally ignores forested areas, which is also a demanding domain in snow hydrological modelling, see Krinner et al. (2018), and references therein. As such, it shares most of the challenges involved in mountain hydrological modelling, which are all related to scaling issues of hydrometeorological conditions in complex topography (Klemes, 1990). In terms of data sources, for past and present meteorological conditions ("nowcast"), in-situ observations can directly be used to drive the models, including, in some cases, snow cover observations. Configurations where NWP model forecasts for past conditions are optimally combined with in-situ observations are referred to as "analysis", consistent with the terminology used in NWP. For future conditions ("forecast"), outputs from NWP models need to be employed and can be adapted in various ways to the geographical model configuration.
In terms of geographical configurations, some model configurations focus on numerical simulations performed "at" observation stations, which are referred to as "station-based" approaches in the following. The alternative is to operate the models on a topography independent from the location of the observation stations, either on a structured grid ("grid-based" approaches) or on the basis of a geometric decomposition of relief ("topographic class-based" approach using mostly elevation, slope angle and aspect descriptors). Fig. 1 provides a graphical synthesis of these approaches.

Currently operational station-based approaches
This section reviews operational implementations of snowpack models using station-based approaches in various countries. This shows Table 1 Overview of driving parameters and conditions for SNOWPACK in operational use.

Type of simulation
Flat field, snow height driven Virtual slope, precipitation driven  Morin, et al. Cold Regions Science and Technology 170 (2020) 102910 commonalities and differences between current implementations. Station-based approaches are the most used, because they provide results, which are considered easiest to implement using in-situ meteorological observations, and to convey to the forecasters as information provided at the location of the observation stations that they are used to.

Switzerland
To support avalanche warning services in their activities, SLF has pioneered operational snowpack modelling using a station-based approach. The operational implementation of SNOWPACK was initiated concurrently with the implementation of a network of automated snow and meteorological stations (nAWS) located in the Swiss Alps. Today SNOWPACK runs simulations at 130 locations whenever new data over at least 1 h is available. SNOWPACK runs in nowcast mode for every S. Morin, et al. Cold Regions Science and Technology 170 (2020) 102910 station from 1 September that year to 31 August the next year. The measured temperature at the snow/soil interface is used as the lower boundary condition (Dirichlet). There is no re-initialization during that period except if driving data are missing for a too long period. In addition, every 3 h and starting with the last output from the nowcast mode, SNOWPACK is driven with the output of a NWP model (COSMO-1, Bellaire et al., 2017). In this case, COSMO data from the closest (Euclidian distance) grid-point are downscaled to the position and elevation of the station. The nAWS are especially designed to work in harsh environment and low power requirement does not allow neither ventilating nor heating of the sensors. Moreover, it is well known how difficult it is to monitor automatically solid precipitation in mountainous environment (Nitu et al., 2018). Accordingly, the input provided by the nAWS include unventilated air temperature and relative humidity, wind speed, reflected short wave radiation, snow surface temperature and snow height, which has two functions. First, it drives the accumulation during snowfall, second allows to correct for errors during subsequent settlement periods. Indeed, those errors lead to wrong estimates of the depth of the next snowfall.

Canada
In Canada, the sparsity of weather stations prompted research into producing "nowcast" simulations using with past NWP model forecasts. Bellaire et al. (2011) and Bellaire and Jamieson (2013) introduced this approach by forcing SNOWPACK with data from the regional scale NWP model from the Meteorological Service of Canada. Weather data is extracted from a set of NWP model grid points and treated as a network of remote weather stations. This approach was operationally implemented in 2010-2011 for Avalanche Canada to assist public avalanche forecasts in remote regions. Simulated snow profiles were produced daily at eight NWP model grid points. Since then, operational upgrades have been implemented on a regular basis including switching from a regional (15 km) to high-resolution (2.5 km) grid in 2013-2014, and adjusting the NWP model data to a standard set of elevations . Since 2016-2017, Avalanche Canada has station-based simulations at 80 grid points throughout southwestern Canada. Simulations are produced at three standard elevations for each grid point (1600, 2000, and 2400 m) by assuming constant lapse rate adjustments for precipitation, air temperature, and humidity. Simulations are updated by linking up 6 h of data from each successive run of the high-resolution Canadian NWP model (High Resolution Deterministic Prediction System, HRDPS). Forecasts for predictive hours 6-12 are used to minimize "spin-up" errors.

Italy
In Italy, the AINEVA avalanche warning services are operationally running the SNOWPACK model since 2004 in the Veneto Region, where it is running operationally in nowcast mode at 25 stations by using AWS data as input. In the Aosta valley, a combination of observed and simulated data (calculated by the NWP model COSMO-2 at 2.2 km resolution) is used for driving the SNOWPACK model in nowcast mode at 6 locations. At the Bormio and Livigno Avalanche Centre (Lombardia region) SNOWPACK simulations are calculated daily for 33 locations by interpolating, using the MeteoIO library, the weather data measured by the existing AWS network located within the area (Monti et al., 2016b). Since winter 2017-2018, both Veneto Region and the Livigno Avalanche Centre run SNOWPACK in forecast mode by using the output of the NWP model COSMO-2 for a total number of 22 locations. In the Marche region (middle of Italy on the Apennines mountain range) SNOWPACK runs for 7 AWS and 7 virtual stations.

Austria
In Austria, the national weather service ZAMG, which operates two of the seven regional avalanche warning services in Austria, recently developed a SNOWPACK-based model chain, which merges station-based "nowcasts" with NWP-based "forecasts" at selected core-stations in mountainous regions (Gobiet et al., 2016). ZAMG currently operates a ensemble-approach by driving the SNOWPACK model with various NWP models at different geographical resolution, in order to quantify the reliability of the forecasts for operational use, although no elevation-dependent adjustment is considered at the moment. In addition, NWP post-processing methods are investigated to improve the representativity of the NWP output for the location of the station.
Since winter 2016-2017 this system is operationally used at the Styrian avalanche service (http://www.lawine-steiermark.at/), with a 6-hour update frequency. The forecasts are driven by five different NWP models and in addition by five different members of the ALADIN-LAEF (Aire Limitée Adaptation Dynamique Développement InterNational -Limited-Area Ensemble Forecasting) NWP ensemble (Wang et al., 2011).

Norway
The Norwegian Avalanche Warning Service has implemented the Crocus snowpack model (within the SURFEX platform) at 82 automatic weather stations with hourly forcing from 2011 to 2016. Since the automatic weather stations lack a standard set of sensors, station data was complemented by forcing data provided by the operational NWP at the Meteorological Service of Norway, e.g. at a station measuring temperature, humidity, and wind speed, precipitation and radiation were taken from the closest grid cell of the NWP. The AROME-MetCoop NWP ) has a 2.5 km grid spacing and is updated four times a day. The 06:00 UTC run was used to drive Crocus once a day. The model was initialized without snow on the 1st of September each year and ran throughout the season without correction until the 31st of May the following year.
This setup was discontinued since the 2017/18 season. Reasons for this were unreliable results and lack of personnel to maintain the model chain. The combination of automatic weather stations with a broad variation of sensors and NWP data as input to Crocus made it hard to monitor and validate the performance of the snowpack model. Model output was often far from observations. The project was therefore put on hold until better model input could be guaranteed.

France
The SAFRAN system described in Section 3.2.1 makes it also possible to extract analysis and forecast time series at arbitrary elevations such as those of observation stations. However it is only since 2015 that such a model configuration has operationally been used for, as of 2019, 665 manual and automated stations in the French mountain regions performing snow height measurements, accounting for the local topography of stations including solar masks from surrounding topography. The latter comes either from direct field measurements, or was computed from a digital elevation model. This configuration operates in both analysis and forecast mode, similar to the the topographic class implementation described below. It is mainly dedicated to a real time evaluation of the model skill.

Currently operational topographic classes approach
Topographic-classes approaches, also termed as "lumped" or semidistributed (Klemes, 1990;Fiddes and Gruber, 2012), are based on a decomposition of terrain features (mostly topographic) designed to sample the most important aspects of land surface heterogeneity without impairing the computational efficiency of the model, especially in data scarce areas where the spatial resolution of available data (be it observations or atmospheric models output) is significantly coarser than the typical length-scales of the local topography.

France
In France, in order to cover a wide range of elevations within the French avalanche forecasting areas (referred to as "massifs"), despite the lack of high-elevation meteorological stations capable of measuring all the data needed to run Crocus operationally, the Météo-France Snow Research Center (CEN) has developed the Système d'Analyse Fournissant des Renseignements Atmosphériques pour la Nivologie (SAFRAN, which stands for "Analysis system providing atmospheric information relevant to snow on the ground", Durand et al., 1993). The key hypothesis behind SAFRAN is the homogeneity of massifs (typically 1000 km 2 ) within which flat field meteorological conditions are assumed to depend only on elevation.
SAFRAN was initially developed in analysis mode (hence its name), and was only later extended for forecasting purpose (Durand et al., 1998). SAFRAN uses as many as possible real-time ground observations for the past 24 h period between 6 h00 UTC D-1 (D-4 from 2019 on) to 6 h00 UTC D (including temperature, relative humidity, wind-speed, precipitation amount, fresh snow thickness and density), together with remotely-sensed cloudiness maps and available radiosondes observations (nowadays playing a minor role). The meteorological data for temperature, relative humidity and wind-speed stems at 3 h time resolution from one representative atmospheric vertical profile output extracted at a coarse horizontal resolution from a given NWP model (ARPEGE) forecast for the same 24 h period for each massif. Forecasted values of temperature, relative humidity, wind speed and precipitation from the NWP model are combined with observations using optimal interpolation and more advanced variational assimilation techniques to produce hourly time series of analyzed meteorological data needed by Crocus. Output data are provided by steps of 300 m elevation (900 m, 1200 m, etc.). For each elevation level, virtual slopes can be generated through modifications of the incoming direct shortwave radiation. The operational configuration of SAFRAN uses flat field, and 20°and 40°s lope angle and 8 slope aspects (N, NE, E, SE, S, SW, S, NW), totalling 17 orientation configurations for each elevation level. In forecast mode, SAFRAN operates in exactly the same way but uses only NWP model output for the periods from D 6 h to D+1 6h and D+1 6h to D+2 6h. Starting in 2019, SAFRAN forecast are provided up to D+4 6 h using the ensemble approach developed by Vernay et al. (2015). This ensures that the coming day is entirely covered by the SAFRAN forecast when the avalanche hazard forecast is prepared on D, valid from D 16h00 to D+1 16h00 approximately.
The system runs for 23 massifs (regions) in the French Alps, 2 in Corsica and 23 in the Pyrenees (11 on the French side, 11 in Spain and 1 for Andorra), and additional forecasting areas in 8 lower elevation forecasting areas in Jura, Vosges and Massif Central. Crocus model runs are performed corresponding to all analysis and forecast configurations of SAFRAN described above.

Canada
In Canada, a topographic class approach has been operational since 2014-2015. The approach is conceptually similar to SAFRAN, which uses data from near surface and free atmosphere elevation levels, but uses purely NWP model surface field forecast data to produce meteorological data for different topography classes. The current implementation divides southwestern Canada into 34 massifs, and within each massif NWP grid points from the HRDPS model (2.5 km) are assigned to below treeline, treeline, and if applicable, alpine vegetation bands based on the elevation of the model grid point. A single SNO-WPACK simulation is produced for each vegetation band by averaging the NWP model data from the grid points within the massif (Horton and Jamieson, 2016). Meteorological data is taken from the same model runs and hours as the Canadian station-based approach.

Currently operational gridded approaches
This section reviews currently operational gridded approaches, which all rely on direct use of NWP model output, with various levels of post-processing and downscaling. Such approaches are the most recently (and actively) developed at present.

Canada
A gridded approach is implemented in Canada with a focus on mapping areas with critical buried surface hoar layers. The model is based on a simplified surface energy balance model that predicts surface hoar formation (Horton et al., 2015, see Section 3.4). The surface hoar formation diagnostic is computed for over 25,000 grid points from the native grid of the HRDPS model (2.5 km grid spacing). The nowcast produced daily with meteorological data from the same model runs and hours as the Canadian station-based approach. This approach was operationally implemented in 2013-2014, and is available to Canadian Avalanche Association members as a suite of maps that show the modelled size of surface hoar for widespread surface hoar layers. Since 2017-2018, the complete SNOWPACK model has been run at all HRDPS grid points in selected study areas (covering over 3000 grid points).

Norway
The main model work-horse for the avalanche warning service in Norway is the seNorge model (Saloranta, 2012(Saloranta, , 2014. The seNorge snow model operates with 1 km grid spacing. It uses gridded observations of daily temperature and precipitation (Lussana et al., 2018) as its forcing for recent days and re-gridded temperature and precipitation from the NWP Arome-MetCoOp  for the next two days and ECMWF for the coming seven days. The seNorge model provides a relatively simple, not very data-demanding, yet process-based method to construct snow maps of high spatio-temporal resolution (1 km grid spacing). It is an especially well suited alternative for operational snow mapping in regions with rugged topography and large spatio-temporal variability in snow conditions, as is the case in the mountainous Norway (Saloranta, 2012). Since it is a single layer model it cannot provide any information related to snowpack stability, such as the existence of weak layers. seNorge and its gridded forcing data are being improved constantly and will soon be available with a temporal resolution of 3 h instead of 24 h.

Austria
The physically-based and spatially distributed snow cover model SNOWGRID (Olefs et al., 2013) is run operationally at ZAMG on a model domain (roughly 400 × 700 km) currently focused on the eastern European Alps and outputs snow data in a very high spatial resolution of 100 m (7001 × 4001 = 28 millions grid points) with a time resolution of 15 min in near real-time and a 72-H forecast.
SNOWGRID is driven with two kinds of gridded meteorological input data depending on its mode of operation: (1) The analysis mode provides output in near real time and uses operational data from INCA (Integrated NowCasting and Analysis system; Haiden et al., 2011), (2) The 72-H forecast mode is based on the NWP models ALARO and AROME (Wang et al., 2006;Seity et al., 2011), running operationally at ZAMG and the ECMWF model.
INCA was developed and is operated by the Austrian weather service ZAMG and uses remote sensing and radar data as well as ground observations providing data on a 1 × 1 km grid with a temporal resolution of 15 min (precipitation, cloudiness) and 1 h (all other parameters). SNOWGRID forecasts are driven by a weather-type dependent weighted average of the ECMWF (8 km grid spacing) and AROME (2.5 km grid spacing) models for precipitation input and the AROME model alone for all other input data (bias-corrected air temperature). Currently all meteorological input grids are simply bilinearly interpolated towards the SNOWGRID grid spacing of 100 m.

United States of America
The Colorado Avalanche Information Center (CAIC) implemented an in-house NWP model in 2011 (Snook et al., 2005;Snook, 2016). The Weather Research and Forecast (WRF) model (Skamarock and Klemp, 2008) was configured with a 4 km grid spacing, which at the time was three times smaller than the highest resolution US National Oceanographic and Atmospheric Administration (NOAA) National Weather Service (NWS) forecast model. CAIC recently decreased the model grid spacing to 2 km, hence, providing better representation of mountainous terrain.
Since many areas of the CAIC forecast zones lack regular observations, the Center started an investigation into driving the SNOWPACK model using WRF model weather forecast data with the objective of providing basic snowpack structure information in these remote areas.
The WRF model is initialized four times per day at 00, 06, 12, and 18 UTC using data obtained from the NOAA. Local weather observations (e.g. CAIC weather station and ski area) are also incorporated. Each WRF model run outputs a complete set of weather forecast information (surface and upper atmosphere) at hourly increments to 84 h. In addition, the model is configured to output surface forecasts every 10-minute with the sole purpose of providing input to the SNOWPACK model.
The WRF model uses a "cold-start" initialization meaning that water vapor information is utilized, but cloud and precipitation information is not available at model startup. The model then takes several forecast hours before the model successfully generates clouds and precipitation, which is referred to as the "spin-up" time. To avoid the spin-up period, forecasts between 6 and 12 h are used as input to the SNOWPACK model. The 6-hour periods are concatenated together from each successive model run to provide a time-continuous data stream at 10minute increments to the SNOWPACK model. Since forecasts are available at every model grid point (2 km grid spacing), then it is possible to run SNOWPACK at each of these points. Snowpack profiles are then updated once per day using the latest 24 hours of weather forecast information.

Slovenia
The In Slovenia, daily gridded SURFEX/Crocus simulations are performed. The domain is a regular 1 km grid with 401 × 301 points centered on Slovenia (the same domain as Inca-SI nowcast that is use operationally). Two daily runs are performed, an analysis and forecast cycle, which differ in the data used to construct the forcing files and the length of the simulation. The execution of the model is implemented using ECMWF's ecFlow workflow package (see Fig. 2). The scripts for constructing forcing files and visualizations are written in R.
In the analysis cycle, A 24 h from 04 UTC -24 h to 04 UTC is performed daily at 04 UTC. The data, used to construct the forcing file is obtained mainly from INCA-SI analysis. The analyzed fields available in INCA are: 2 m temperature and relative humidity, 10 m wind speed and direction, 30 min total precipitation, precipitation type (rain, snow, or mix) and short wave radiation. The long wave radiation also needed to run SURFEX is obtained from the Aladin NWP model. The forecast run is performed after the analysis run is complete. The forcing files are constructed from the 00 run of the Aladin NWP model, which operates on a 4 km grid. The relevant fields are transformed to the 1 km grid by a simple bilinear interpolation. The forecast range of the NWP model is 72 h, which yields a 68 h forecast of snow conditions.

Post-processing of snowpack model outputs
Raw 1D model output consists of the vertical profile of simulated physical properties of snow (one value for each property and each numerical layer), along with diagnostics of the energy and mass fluxes at the interfaces. Such information needs to be post-processed into process-relevant information which can be used for avalanche hazard forecast activities. Below we describe the post-processing steps currently operationally used downstream the snowpack modelling chains, mostly at the point scale but also in a few cases at the regional scale (see Section 2.6). The selection of post-processed variables made available for the forecasters and their visualization are given in Section 6.

Simple diagnostics
Total snow height and snow water equivalent (SWE) are a typical diagnostic of physically-based snowpack models. Snow surface temperature and bottom liquid water flux also inform on the thermal state of the snowpack.
Snowpack models described above can be used to derive the total thickness and mass of snow deposited for a given time period (e.g., 1 day, 3 days etc.). This allows for computing a recent snow amount (depth of snowfall and its water equivalent) from the model output, which accounts for settling of freshly fallen snow and thus can be compared to snow board measurements (Fierz et al., 2009). In addition, multi-layer models make it possible to compute the near surface wet snow thickness (sum of the thicknesses of contiguous layers with a non-zero liquid water content, starting from the top if the uppermost layer is also wet), as well as near surface refrozen thickness (sum of the thicknesses of contiguous layers with a null liquid water content, starting from the top if the uppermost layer is also dry).

Mechanical stability diagnostic
While Crocus only addresses physical transformation processes occurring in the snowpack, diagnosed mechanical properties and snow stability have been handled by the MEPRA model (Giraud, 1992). MEPRA computes the penetration resistance and shear strength of each layer depending on density, snow type descriptors and thermal state. Furthermore, depending on the vertical profile of physical properties, MEPRA assesses the location of potential weak layers and computes stability criteria for natural releases and skier triggered avalanches for each profile simulated by Crocus on 40°slopes. The computation of the avalanche hazard index depends not only on the value of the current instability criteria but also on their variations in time. For example, if the mechanical stability undergoes an increase, the diagnosed hazard level will be lower than the exact same current situation but with a decreasing trend of the mechanical stability.
SNOWPACK features post-processing algorithms for directly providing snow stability information. Schweizer et al. (2006) implemented the snow layer hardness parameterization by performing a multivariate statistical regression analysis for the snow hardness index for each snow type using snow density and grain type as independent variables. Monti et al. (2014) improved this parameterization by implementing the snow hardness as a discrete variable. The dry snow stability evaluation within SNOWPACK profiles is performed in three independent stages. First of all, potential weak layers are detected using a threshold sum approach inspired by manual snow profile interpretation systems (Schweizer et al., 2008;Monti and Schweizer, 2013). Once the potential weak layers are detected, their strength regarding the fracture initiation is evaluated by using a multi-layered skier-induced stress approach (Monti et al., 2016a), which is a refined version of the classic skier stability index (Föhn, 1987;Jamieson and Johnston, 1998) (see Fig. 5).
A third index provides information on fracture propagation in weak layers, an essential process for avalanche release (Gaume et al., 2017). This index was only recently implemented and still requires a more extensive verification. Mitterer et al. (2013) introduced an index describing the wetting state of the snowpack, which is operationally used in nowcast and forecast mode to assess the advance of the snowpack wetting based on SNOWPACK model runs (since 2011/2012 at SLF). This index is defined as:

Wet snow diagnostics
where Θ w, v represents the average liquid water content of the full snowpack. Even though this approach strongly simplifies the wetting of the snowpack, it is easy to interpret by forecasters. Strong increases in the value of the median LWC index often correspond with the onset of wet snow avalanche activity (Mitterer et al., 2013;Bellaire et al., 2017).

Surface hoar diagnostic
Horton et al. (2015) proposed a method to predict surface hoar formation. It is based on a simplified surface energy balance model, which feeds the surface hoar formation routine from SNOWPACK (Lehning et al., 2002b). This allows for the computation of the size of surface hoar at the time of burial. The accumulated precipitation since the time of burial is also calculated at each simulation point to estimate the load on buried surface hoar layers.

Drifting snow diagnostic
SYTRON (Durand et al., 2001;Vionnet et al., 2018) and the drifting snow module of SNOWPACK (Lehning and Fierz, 2008) provide indices of drifting snow (occurrence and intensity) in a topographic class and station-based approaches, respectively.

Post-processing of snowpack model outputs at the regional scale
Post-processed information at the point scale can be aggregated at a large scale in order to provide synthetic information relevant to a larger area. For example, in Switzerland the point-scale LWC index is aggregated as a function of elevation bands, thereby providing an regional wetness index, based on the median of individual LWC index values, reducing the potential bias or error of individual stations (see Fig. 10, right).
In France, the mechanical stability indices for natural avalanche release calculated by MEPRA for each simulated snowpack profile is aggregated at the massif-level, in an attempt to provide a massif-scale natural avalanche hazard level (Vernay et al., 2015). This index ranges from 0 (lowest level) to 8 (highest level), and is built as a combination of MEPRA natural hazard indices at different elevations from 1500 to 3000 m for 40°slopes and 8 aspects. It was designed to be comparable to a massif scale index of observed avalanche activity (Giraud et al., 1987). An integrated index for triggered avalanche hazard has been developed but never thoroughly evaluated and made available to the forecasters. In theory, such indices, computed at the regional scale, which in the case of France correspond to the avalanche hazard rating zones (Techel et al., 2018), bridge the gap to the regional hazard level assessment. However, they do not account for the simulations including wind-induced snow transport.

Technical challenges posed by operational use of snowpack models
It is assumed here that model developers are functionally separated from their users. This implies that automated model chains have to be used, targeting routine operations regardless of the date and the day of the week during the course of the snow season. Model developers may not be available in case technical questions (including potential shutdowns) arise during the course of operations, e.g. during holiday or week-end periods, or very early in the morning.
In contrast to research developments which usually do not require real time operations of models, running snowpack models in support of operational forecasting activities requires that appropriate measures are taken to ensure continuous operations of models, in two respects. Firstly, if the model results are part of the information routinely used by the forecasters in their daily duties, then the model needs to be able to operate without interruption and continuously provide data to be visualized and analyzed. Adequate measures are to be taken to fix potential discontinuations of model chains, be it for intrinsic model errors or for external reasons such as diverse causes of failure of the information technology system hosting the model chain (operating system stability, network connection and power reliability, delay in input data access, failure of database requests etc.). Breaks in the data flow provided to the forecasters can be prejudicial to their daily activities. Secondly and perhaps more importantly, because of the long term (weeks to months) memory of the snowpack and the fact that incremental changes of its status occur throughout the snow season, breaks in the data stream to the snowpack model itself could have adverse consequences not only for the period of time considered by the data gap but for the rest of the snow season unless appropriate mitigation is implemented. For example, if observation data directly used as input for the snowpack model are missing or worse, incorrect, this may halt or at least significantly impair the simulation process for the rest of the snow season.
Three main approaches are taken to handle the first category of issues, i.e. those related to supplying continuously simulation results. First of all, all the codes used should be written in the most possible robust way, and tested thoroughly to avoid failures due to programming errors. Second, a stable IT framework should be used for the implementation of the operational model chains. As an example, Fig. 2 shows the workflow used for the implementation of Crocus in Slovenia, which uses the European Center for Medium Range Forecasts (ECMWF) operational tool ecFlow. Third, continuous real-time monitoring of model tasks should ideally be performed, in order to identify model errors and re-launch the model chains in case of externally-driven failure (e.g. due to temporary data unavailability). This implies that the model chain technical description is sufficiently detailed so that personnel in charge of monitoring the operations of model chains are able to perform their duties without direct support from the model developers.
Fixing the second category of issues (i.e., those related to the absence of forcing data and how this impacts the continuity of the snowpack simulation during the course of the snow season) is highly model dependent. As a general rule, robust and redundant observation data streams help reducing the number of data gaps. Pragmatic automated gap-filling techniques, e.g. using data from neighbouring stations, can lead to enhanced resilience of the model chains to gaps in the observation data streams, at the expense of the accuracy of the model results. Finding the right balance between accuracy and continuity of the simulations can be challenging, but is an inherent component of the design of the model chain.
For example, Météo-France snowpack model chains are constantly (24/7) monitored by the technical personnel also in charge of monitoring NWP model chains, and can act upon the unavailability of expected products from the model chains, or the issuance of error diagnostics. Ideally, error diagnostics should be handled by the model chain itself and their status should trigger further automated trials of the failed tasks. In the Météo-France chain of models, the main source of failure is due to the temporary unavailability of either NWP forecasts or observation datasets, which are used by SAFRAN, and failure to access the database containing the initial state of the snowpack (from previous day simulations) necessary to restart the model the following day. Major modelling errors (i.e. causing a stop of the model chain, which can only be fixed through a modification of one or several of the programs, and may lead to the absence of critical data for the continuity of the snowpack simulation through the winter season) have not been encountered in several decades of operations. This strong robustness of the model chain however has a downside: it has occurred several times that a minor error (i.e. not leading to model failure) systematically affected the model chain but was only noticed several weeks after their effect, because of the difficulty for model users to disentangle true programming errors or model failure from model errors arising from the imperfectness of numerical models.
Ways to identify minor technical failures of model chains may develop twofold, based on experience gathered in operational NWP. First of all, model chains should be designed to automatically issue operational logs highlighting any unforeseen behaviour of the programs (e.g., monitoring of the quality/quantity of observation or NWP data used by the model). Secondly, quantitative performance metrics should be computed operationally and monitored, in order to rapidly address cases when the performance strongly deviates from the norm.
Technical failures of the modelling system are reduced by an automated logging and alert system in the whole model chain including the post processing of products and by constantly computing model performance and plotting time-series of model vs. observed snow height values at all main evaluation sites in a visual environment where archived model data can also be accessed (see Fig. 3). Operational NWP organizations typically host teams specializing in monitoring such information, and providing feedback to model developers and IT departments as soon as abnormal situations occur. Even if avalanche forecasting centers may not afford such investments, good practices from operational NWP organizations could still serve as inspiration regarding the design of model chains and their monitoring component (Commission for Basic Systems, 2002).

Information quality of snowpack model output
Previous sections have addressed existing model chains and their operating configurations. We now turn on the information quality of the model output, taking the viewpoint of the avalanche forecasters. We use the information quality framework of Bovee et al. (2003) to describe the usefulness attributes to the snowpack model outputs. Bovee et al. (2003) state that "to determine the quality of information we must (1) get information that we might find useful (accessibility), (2) understand it and find meaning in it (interpretability), (3) find it applicable to our domain and purpose of interest in a given context (relevance), and (4) believe it to be free from defects (integrity). We would dismiss or discount information that meets our criteria for all but one of any of the aforementioned attributes, each of which may be more than just a binary value." Below we outline how these four attributes can be substantiated in the case of snowpack model output for avalanche forecasting. This brings together knowledge from the developers of model chains and from the users community. General feedback on the use of operational forecasters has been gathered from the organizations which operate snowpack models in support of avalanche hazard prediction. A detailed, panel-based evaluation of their use of this resource is beyond the scope of this article, nevertheless recurring topics were identified and are outlined below where relevant.

Accessibility
Accessing snowpack model outputs is obviously required for their use by the forecasters. This can be decomposed in two aspects, time and simplicity.
• Time. Avalanche hazard forecasters work in a time constrained environment. Information from various sources, sometimes spanning large geographical areas, need to be processed and analyzed to produce a nowcast or forecast of avalanche hazard level for various rating regions. In many cases, the information from snowpack models is seen as additional to the more traditional sources of informations (observations and weather forecasts), and under time constraints the use of snowpack model outputs tends to be considered only if time allows. Full integration of model output analysis into an existing workstation helps reducing this trend, but does not fully alleviate it.
• Simplicity. Uptake of snowpack model output by operational avalanche hazard forecasters was enhanced in case where a limited number of products, most often directly analogous to otherwise observed information (e.g. stratigraphy, height of new snow etc.). In France, despite the provision of regional-level integrated avalanche hazard level predictions, and a GUI which favors progressive deepening of the analysis level from the massif scale, to the elevation/ aspect scale, finally to the individual stratigraphic profile as the last stage of the complexity level, forecasters' most common feeling with the model output is to be overwhelmed by the quantity of information made available to them. Although it could be useful for advanced users, excessive amount of information tends to reduce the overall accessibility of the model output for the general population of forecasters. S. Morin, et al. Cold Regions Science and Technology 170 (2020)

Interpretability
Model output needs to be interpretable by the forecasters to be useful. Interpretability is maximized when model output corresponds to an otherwise observed quantity and is visualized in a manner comparable to the visualization of observations. Training and effective visualizations are enablers for the interpretability of the model output.

Relevance
One key sub-attribute of relevance, according to Bovee et al. (2003), is the so-called datedness of the information. As long as the data exists and is accessible, a time tag is generally provided along with the information, so that this attribute, although critical, is generally not a concern for operational avalanche forecasting. More interestingly and critical, it is often questioned, whether snowpack model output provide forecasters with information, which is not already available directly or indirectly through other data sources (observations, weather forecast etc.). Addressing this question quantitatively is difficult, and would require conducting experiments with several parallel teams of forecasters over the same avalanche situation to assess the added-value of each piece of information.

Integrity
In the information quality framework of Bovee et al. (2003), information integrity is composed of accuracy, completeness, consistency and existence. Completeness, existence and datedness can be addressed simultaneously to the accessibility attribute. Accuracy and consistency are key attributes, which generally hamper the overall credibility, hence   The markers within the snow cover represent the depth of potential weak layers detected by using the relative threshold sum Approach (Monti et al., 2014). The colors of the markers are related to the values of the multilayered skier stability index (SK ml ) (red = SK ml ≤ 0.5; orange = 0.5 < SK ml ≤ 0.95; yellow = 0.95 < SK ml ≤ 1.05; light green = 1.05 < SK ml ≤ 1,5; green = SK ml > 1.5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) S. Morin, et al. Cold Regions Science and Technology 170 (2020) 102910 usefulness, of the model outputs. Accuracy, in this framework, "refers to information being true or error free with respect to some known, designated or measured value". Operational use of numerical models has to cope with model errors, in particular the significant misses of the models. These are unfortunately unavoidable, and can greatly hamper the credence given to the models (Pappenberger et al., 2011). The occurrences of significant perceived errors in precipitation amounts, for example, and the resulting snow height, is a typical case where avalanche hazard forecasters tend to reduce their trust in the model outputs, even in cases where useful information can still be drawn from the model results. Verification of snowpack model output is challenging for a variety of reasons, which can all be traced to the various sources of uncertainty which affect every component of the model chains. Uncertainties are associated to intrinsic observation and forecast errors at the point scale (applicable both for meteorological and snow data), but also the spatial representativeness of data used for the evaluation of model results, be it at the point scale (Grunewald et al., 2013) or remotely-sensed. In fact, there has been little systematic undertaking of the verification of the predictive capability of physically-based snowpack models and their associated post-processing routines in terms of avalanche hazard forecast due to the absence of objective field measurements of avalanche Fig. 7. Numerical simulations of snowpack profiles spanning several elevation levels for a given massif, aspect and slope, from the SAFRAN-Crocus-MEPRA model chain used in France (example from Haute-Bigorre, Pyrénées, West-facing 40°slope). Stratigraphy data include temperature (blue), density (green), mechanical resistance (brown), liquid water content (red) and snow type (colour coding according to Fierz et al. (2009)). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) hazard at the regional scale (Giraud et al., 1987;Bellaire et al., 2017). This also applies to snowpack stability, which has lacked hitherto objective and reproducible field measurements (Reuter et al., 2015). One step less in the model chain, regarding snow stratigraphy, comparisons between observed and simulated stratigraphy data have only been carried out for a few selected cases (Brun et al., 1992;Durand et al., 1999;Lehning et al., 2002a) and generally lacked an objective framework for quantitative comparisons between observed and simulated data in a context with discrepancies in layering between observed and simulated data are the norm Fierz and Lehning, 2004;Hagenmuller and Pilloix, 2016;Hagenmuller et al., 2018). What has mostly been carried out consists in comparing simulated and observed snow height and, when available, snow water equivalent Lafaysse et al., 2013;Schmucki et al., 2014;Bellaire et al., 2011). However, these general evaluation studies are often performed irregularly in time (often when a new version of the model chain is released), are not necessarily communicated clearly to the users, so that they do not contribute to increasing the credibility of  The lines show the LWC index from the retrospective SNOWPACK simulation for a south-facing slope with a slope angle of 38 degrees. The blue line represents the nowcast (i.e. SNOWPACK is driven by meteorological observations) and the coloured lines represent the forecasts from different NWP models (i.e. SNOWPACK is driven by an ensemble of NWP models). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) model results. As described in section 5, only continuous performance monitoring makes it possible for users to infer, in real time, the performance level of the model chain under ongoing meteorological and snow conditions -given that the performance level can vary in time, depending on snow conditions, for example. The matter is further complicated by the fact that in-situ observations are affected by horizontal variability, questioning their representativeness, so that metrics computed using simulated and observed quantities can be challenged on both sides (simulations and observations) by users.
Consistency of model output is intrinsically perfect, because all model outputs are generated from the same model run and are therefore produced from an inherently consistent computational process. Perceived inconsistencies can either stem from incomplete access to model output (thereby impairing the user's ability to identify the reasons behind seemingly spurious evolution of several model outputs), or from the user's inability to understand such changes because of

Examples of graphical representations of model outputs provided to forecasters
In this section we illustrate typical graphical representations of model outputs provided to forecasters.

Snowpack profiles
Snowpack profiles are probably the most iconic outputs of physically-based snowpack models. They consist of the stratigraphic representation of the simulated layering of snow properties (snow type, density, temperature, liquid water content, grain size, shear resistance) and sometimes diagnosed critical layers (weak layers). Snowpack profile information can be displayed in several ways. Fierz et al. (2009) provides a series of examples of such representations for observed stratigraphic profiles. Below we illustrate examples from simulated profiles. Fig. 4 illustrates the time evolution of snow type, simulated using SNOWPACK driven by in-situ observations from the Weissfluhjoch field site near Davos. Fig. 5 illustrates the visualization of indicators within the snowpack (SNOWPACK implementation in Italy), directly identifying weak layers and their corresponding mechanical stability. While the previous examples apply regardless of the geometrical configuration used for the model chain, snowpack profile information can be specifically plotted using a station-based approach (Fig. 6) or a topographical class-based approach (Fig. 7). Both representations provide direct vizualisation of simplified snowpack stratigraphy information depending on geographical location (Fig. 6) or for a given massif, elevation and slope angle, depending on the slope aspect (Fig. 7, from the SAFRAN-Crocus model chain using the Synopsis visualization system, see Coléou et al., 2016). Fig. 8 shows multiple simulated profiles from a gridded approach by stacking the stratigraphy of a single snow property (snow type) side-by-side (Horton et al., 2018). This alternative to traditional snow profiles could allow general overviews of the main trends in large batch simulations.

Scalar information
Scalar information is both input to the snowpack model (meteorological forcing data for past, present and future conditions) and output (post-processed diagnostics, see Section 4). Regardless of the geometrical configuration of the model chains, scalar data can typically be visualized as time series for one or several simulation points, allowing direct comparison of observed and simulated information, for past, present and future conditions.

Time series
We illustrate here the plotting of scalar information in the case of the LWC index introduced by Mitterer et al. (2013) (see Section 4.3) as implemented by the ZAMG for the avalanche service of Styria (see Section 3.1.4). Fig. 9 shows an example of increasing liquid water The LWC index is presented in two ways at SLF, see Fig. 10, either as an elevation plot (Fig. 10, left), allowing the elevational estimate of the wetting front, or a plot, showing the temporal evolution at three different elevations (Fig. 10, right). The latter plot uses an ensemble-type wetness index (median LWC index ), reducing the potential bias or error of individual stations.

Pie plots/aspect-elevation roses
In the topographic-classes based approach, scalar information can be represented as pie plots (also referred to as aspect-elevation rose). This is well suited to identify variations with elevation, slope angle and aspect. Fig. 11 shows an example of a pie plot (aspect-elevation rose) for the natural and accidental hazard level diagnosed by MEPRA from SAFRAN/Crocus simulations, along with a time series and vertical profile for one topographic configuration. Fig. 12 shows an example of visualization of data from the snowdrift modelling system SYTRON implemented within the topographicclass geometrical approach of the SAFRAN-based modelling chain in France . Although much more detailed information could have been provided, for example snowpack profiles corresponding to each simulation points, based on interactions with forecasters a relatively simple visualization chart was designed and implemented, focusing only on the intensity of the wind-driven erosion and deposition patterns.

Maps
Station-based approaches make it possible to draw maps where individual scalar values at one point in time are represented. Fig. 13 shows an example where individual values of simulated and observed height of new snow are represented, along with coloured overlay obtained by spatially interpolating individual values. Such a representation, however, does not represent topography explicitly.
Grid-based approaches are ideally suited to produce continuous maps superimposed on topographic information. In Canada, model data are displayed as map overlays that are viewed on the same basemap as other geospatial avalanche data including avalanche observations and regional hazard ratings. Examples of map overlays include colour shading to represent regions with modelled surface hoar and contour lines to represent modelled snow water equivalent. Fig. 14, from Horton et al. (2014), shows an example of a map of surface hoar formation together with an estimate of the load of new snow upon it. However, it should be kept in mind that at high resolution, the representation of variables strongly influenced by the topography (e.g., snow height) can be very noisy and difficult to interpret from a regional perspective. Fig. 15 shows an example of a map displaying the penetration of a ram sonde tube, indicative of the mechanical strength of surface snow, for the entire geographical domain of Slovenia. Fig. 16 shows an example of series of maps of the regional averaged avalanche danger level diagnosed by MEPRA for various dates in the French Alps, using the Synopsis software (Coléou et al., 2016).

Graphical User Interface and data overlaying
All of the representations introduced above are only useful for the forecasting activities if navigation across information types (observations, model results, meteorological forecast) and representation methods (profiles, time series, pie-plots, maps etc.) is as intuitive as possible. This highlights the need for efficient Graphical User Interfaces (GUIs) taking into account both the users needs and the technical capabilities of the model chains used. Powerful and integrated GUIs allow for in-depth analysis of a large variety of information through efficient overlaying of different data sources onto the same graphical representation. Fig. 13 is a good example showing the integration of observed and simulated height of new snow from various data sources. Fig. 17 shows an example of SAFRAN/Crocus data at the massif scale (here cumulated precipitation at a given elevation), overlaid on the output of the NWP model ARPEGE. This allows the user to infer the potential discrepancies between various meterological sources at different space and time resolution, one of whom being used to drive the snowpack model (Coléou et al., 2016).

Future challenges
Improving the scientific and technical methods for snowpack model chains is critical to improve the quality and quantity of information made available to the forecasters as well as feeding products directly delivered to end-users. Work is also critically needed in the field of data visualization, and the ergonomics of how the data is presented to forecasters. Last, perspectives are discussed for more integrated ways forward, in terms of transnational implementation of models and tools.

Re-assessing how model results are made available to the forecasters
The concerns highlighted above require in-depth re-assessments on how the model results are made available to the forecasters, and how they are used. This issue is not only challenging for avalanche hazard prediction, but more generally to all decision-making processes involving interactions between human forecasters and predictive models (Pagano et al., 2014). While in some operational forecasting domains, e.g. meteorology, the role of human forecasters is constantly reducing in favor of NWP outputs (Sills, 2009), in avalanche hazard forecasting the issue at hand is clearly how to make the most of information available to the forecasters, who in some case rely only marginally, if at all, on numerical simulations in support of their activities. The relevance of providing snowpack model output at an insufficient level of performance is also questionable. Indeed, providing unreliable information can disturb forecasters activities and reduce their possible confidence in improved systems in the future, although feedback from operational services is generally useful. The possibility to automatically produce avalanche warning information remains a long term, possibly elusive challenge (Floyer et al., 2016).
Improvements of the situation may be favoured by the following directions: • Direct research towards adding value to forecasts, through better design and visualization of the post-processed outputs at the point and regional scale.
• Better ways to communicate uncertainty. To add value in a human forecasting chain, models need to reduce the forecaster's perceived uncertainty about their knowledge of snow conditions. However snow models have several levels of uncertainty (meteorological inputs, snow science models, spatial representativeness) and most products are poor at communicating that.
• Better integration of the model outputs with existing data management and visualization platforms would make the products easier to use and more valuable to forecasters. Building systems where model output can be directly compared with field observations is an promising way to verify the model and for users to better understand the performance and value of the model.
• Better combination of snowpack models with other models such as statistical models, additional weather products, and risk-based models, such as the conceptual model of avalanche hazard (Floyer et al., 2016;Statham et al., 2010Statham et al., , 2018 or the auxiliary avalanche warning matrices (Müller et al., 2016).

Improvements on the atmosphere and snow sciences side
Main research items, which are likely to lead to improved snow model chains, can be summarized as follows: • Snow physics and snow modelling. This concerns the improvement of the modelling of intrinsic snow processes, such as liquid water dynamics, snow metamorphism, the impact of light absorbing impurities etc. Advances in snow mechanics could translate into renewed methods to estimate snow stability from simulated snow profiles (Reuter et al., 2015).
• Verification. Improved measurements methods, and recently developed method making it possible to pair observations and simulations despite mismatches in layering (Hagenmuller and Pilloix, 2016;Hagenmuller et al., 2018), open the possibility to more quantitatively evaluate the output of snowpack models and guide most required improvements. Such verification developments may also benefit from remote sensing data, including detection of avalanche activity (Eckerstorfer et al., 2017).
• Data assimilation and ensemble analysis and forecast. Deviations of model output with in-situ and remotely-sensed observations are a critical issue hampering operational use of the models by the forecasters. Data assimilation is a promising, yet challenging way forward in this respect. Recent developments in data assimilation rely on ensemble-based methods (Lafaysse et al., 2017;Cluzet et al., 2020), which, in addition to estimating forcing and model errors quantitatively within the assimilation system, open new possibilities for ensemble analyses and forecasts of avalanche hazard (Vernay et al., 2015).
• Higher spatial resolution. Numerous snow processes, e.g. blowing snow and preferential deposition, operate at the scale of a few meters and in 3D, and thus require an explicit representation of topography and associated processes. Operational blowing snow diagnostics are currently available for virtual crests at the station scale (Lehning and Fierz, 2008) or per elevation band within each massif (Durand et al., 2001;Vionnet et al., 2018). This conceptual representation of the topography fails at reproducing the local and complex features of wind field in alpine terrain. To overcome this limitation, distributed blowing snow models running on high-resolution grid (10-250 m) are available (e.g. Liston et al., 2007;Vionnet et al., 2014). So far, they have only be used for research purpose. Their operational implementation over large mountainous regions to produce daily analysis and forecast is still a challenge due to the high computational times of these systems and to the difficulty to obtain reliable high-resolution precipitation and wind fields. On this point, improvements could be obtained from the joint use of sub-kilometer high-resolution NWP systems in mountainous terrain (Vionnet et al., 2015) and statistical downscaling methods relying on terrain-based parameters Winstral et al., 2017). NWP systems do not only provide information on mean wind speed but includes also parameterization of wind gust (e.g. Schreur and Geertsema, 2008) to derive maximum wind speed. Such information could be used to determine more accurately the occurrence and intensity of blowing snow under S. Morin, et al. Cold Regions Science and Technology 170 (2020) 102910 gusty atmospheric conditions (Naaim-Bouvet et al., 2011). Overall, it is expected that progress in this area can benefit from interactions with developments in the field of mountain hydrology (Marsh et al., 2018).
• Post-processing of model output. Model output statistics are routinely used in NWP and hydrological prediction (Vannitsem et al., 2018), yet seldom used for avalanche hazard forecast, although similar issues need to be addressed. Furthermore, the skyrocketing development of data science and artificial intelligence applied to a broad range of practical issues makes it possible to envision applications to avalanche hazard forecasting (see e.g. Nousu et al., 2019), thereby providing added-value products beyond raw snowpack model outputs, with potential for directly feeding the production of avalanche bulletins.

Technical and organizational perspectives
In contrast to the situation prevailing a few years ago, most snowpack models used by operational services are now open source, community-based software with version control and increasing verification routines. This increases the overall robustness of the models, although the costs associated to their maintenance and development is much higher on their host institutions. Moving further into this direction is a clear way forward, following previous examples from the climate, meteorological or hydrological communities.
Beyond developing and maintaining the codes, it may be worth exploring an increased level of pooling of computational resources, especially for organizations operating in neighbouring geographical domains. Currently implemented station-based or topographic classbased approaches are generally limited to national or regional boundaries in their use. The development of grid-based approaches based on NWP data, remote sensing data and high performance, ensemble based forecast and data assimilation methods, may develop more efficiently in a cross-boundary context, in close association to high performance computing infrastructures, rather than in multiple local implementations. Indeed, many NWP models are currently implemented beyond national boundaries, see e.g. Fig. 17. This approach would facilitate multi-model ensemble systems (multiple NWP models, multiple snowpack models), analogous to meteorological and hydrological prediction frameworks. Such a long term endeavor remains to be defined, but could benefit from multi-national opportunities such as the Copernicus services, in Europe (e.g., Copernicus Emergency Management Services, or Copernicus Land Monitoring Services), or large scale similar initiatives in North America. Such an approach may even make it possible, in the longer term, to deploy snowpack modelling chains over large geographical areas currently devoid of such systems, such as High Mountains of Asia or South America. Given the strong relationships with mountain hydrological prediction context and goals, enhanced collaboration with operational hydro-meteorological organizations could also lead to significant progress, while mitigating potential IT Fig. 16. Example of series of maps of the regional averaged avalanche danger level diagnosed by MEPRA (Vernay et al., 2015) every 3 h in the French Alps.
S. Morin, et al. Cold Regions Science and Technology 170 (2020) weaknesses of individual small forecasting centers.

Conclusions
Early warning systems for a range of weather-related natural hazards increasingly rely on numerical models utilizing numerical weather prediction (NWP) output. This article summarizes currently implemented approaches for snowpack modelling in support of avalanche hazard forecasting in Europe and North America. It shows the similarity in approaches for implementing snowpack models for operational purposes, progressively moving away from pure nowcast approaches based on observations, into increased use of NWP output for future predictions. It also highlights a large diversity in the way the information is post-processed and conveyed to forecasters. The information post-processing and visualization component is probably the weakest point of current model chains for operational avalanche forecasting, hampering their use by the forecasters in many operational centers. Innovative approaches, bringing together modellers, forecasters and communication experts, have the potential to make significant improvements in this area. The generalization of modelling approaches directly using numerical weather prediction output, which cover large regional domains, represents an opportunity to bring together modelling groups from neighbouring areas and collaboratively develop the next generation of model chains in a transboundary context, bringing together specific expertises in snowpack modelling, NWP output downscaling and data assimilation, the combination of whom is critical for further progress, and requiring robust and numerically-efficient computing environments for operational implementations. In parallel to these developments, future research should analyze in more details the workflow of operational forecasters, to better assess the positioning of snowpack model information within all streams of data available to the forecasters. Such an assessment could make it possible to derive recommendations which would benefit both the forecasters and the model developers in their respective tasks. worldwide community. CNRM/CEN is part of LabEx OSUG@2020. We thank Marie Dumont (CNRM/CEN) for feedback on the manuscript. This publication has benefitted from funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 730203. Open-access publication fees were covered by Météo-France. S. Horton was financially supported by Avalanche Canada and the NSERC IRC in Avalanche Risk Management at Simon Fraser University. Comments and suggestions from Karl Birkeland and one anonymous reviewer were useful for improving the original manuscript.