S3M 5.1: a distributed cryospheric model with dry and wet snow, data assimilation, glacier mass balance, and debris-driven melt

. By shifting winter precipitation into summer freshet, the cryosphere supports life across the world. The sensitivity of this mechanism to climate and the role played by the cryosphere in the Earth energy budget have motivated the development of a broad spectrum of predictive models. Such models represent seasonal snow and glaciers with various complexities, and generally are not integrated with hydrologic models describing the fate of meltwater through the hydrologic budget. We 5 present S3M v5.1, a spatially explicit and hydrology-oriented cryospheric model that successfully simulates seasonal snow and glacier evolution through time and that can be natively coupled with distributed hydrologic models. Model physics include precipitation-phase partitioning, snow and glacier mass balances, snow rheology and hydraulics, a hybrid temperature-index and radiation-driven melt parametrization, and a data-assimilation protocol. Comparatively novel aspects of S3M are an explicit representation of the spatial patterns of snow liquid-water content, the implementation of the ∆ h parametrization for 10 distributed ice-thickness change, and the inclusion of a distributed debris-driven melt factor. Focusing on its operational implementation in the north-western Italian Alps, we show that S3M provides robust predictions of the snow and glacier mass balances at multiple scales, thus delivering the necessary information to support real-world hydrologic operations. S3M is well suited for both operational ﬂood forecasting


Introduction
The cryosphere is a decisive driver of the Earth system (Barry, 2011;Beniston et al., 2018).Besides altering surface albedo and so concurring to the regulation of global temperature (Flanner et al., 2011), snow and glaciers accumulate winter precipitation and release it during the warm, spring and summer seasons, when demand by societies and ecosystem services is comparatively high (Barnett et al., 2005).This shift in water supply supports water, food, and energy security across climates (Viviroli et al., 2007), with key societal implications (Sturm et al., 2017).For example, snow represents up to 80% of annual water supply in the semi-arid, largely summer-dry western US (Bales et al., 2006;Serreze et al., 1999;Skiles et al., 2018), while 1.4+ billion people in Asia rely on discharge from high-mountain, cryosphere-dominated regions (Immerzeel et al., 2010).Meanwhile, the Andean cryosphere acts as a significant freshwater resource for semi-arid regions of South America (Masiokas et al., 2020), with an estimated contribution of up to 27% to dry-season water supply in La Paz, Bolivia (Soruco et al., 2015).
Seasonality between winter accumulation and summer melt (Barnett et al., 2005), compounded by equally complex but more short-term processes such as rain-on-snow (Rössler et al., 2014), challenges decision makers like water-resources or hydropower managers, who need early and diverse information about snow and glacier-ice amount, distribution, and melt timing to make accurate decisions on water use, allocation, and storage (Georgakakos et al., 2004;Anghileri et al., 2016;Avanzi et al., 2018).This need has catalyzed the development of a variety of models to predict snowmelt-and icemelt-driven discharge (DeWalle and Rango, 2011), to the extent that cryosphere modeling is a dominating topic of both basic and applied contemporary geosciences (Dozier et al., 2016).Applications for cryosphere models are not limited to water-supply and flood forecasting, but include avalanche hazard forecasting (Bartelt and Lehning, 2002;Vionnet et al., 2012), land-surface and weather modeling (Dutra et al., 2010;Wang et al., 2017), and snowmaking (Hanzer et al., 2020).The projected rise in future temperature and aridity (IPCC, 2021) further prioritizes robust predictions of cryospheric water resources, because shrinking glaciers and decreasing snow accumulation may endanger water supply and its predictability (Harrison and Bales, 2016) especially during droughts (Huning and AghaKouchak, 2020).
While detail in process representation may appear as the prime driver of model selection, in hydrologic practice this choice also depends on other four pragmatic factors, which make hydrology-oriented cryospheric models essentially different from those oriented to, e.g., avalanche-hazard forecasting.First, streamflow generation in cold regions involves not only snow and glacier ice, but also precipitation-topography interactions (Blanchet et al., 2009;Mott et al., 2014;Cui et al., 2020), vegetationwater feedback mechanisms (Zheng et al., 2016;Avanzi et al., 2020), and soil-water storage (Bales et al., 2011).Thus, snow and ice hydrology is inherently spatially distributed and multi-scale (Dozier et al., 2016), with the focus being arguably more on distributed than on point predictions (Blöschl, 1999).Second, high-elevation cryospheric regions remain largely ungauged (Rasmussen et al., 2012;Avanzi et al., 2021), meaning that the necessary input data to run complex models is often missing or sparse at best.This condition has favored parsimonious models (Bartolini et al., 2011) and data-assimilation schemes to remedy model deficiencies with independently observed data (Andreadis and Lettenmaier, 2006;Piazzi et al., 2018).Coupled with data sparsity is the third factor, that is, the evidence that simplified and complex models often yield comparable predictive accuracy for bulk processes relevant to the seasonal freshet, such as Snow and Ice Water Equivalent (Huss et al., 2010;Avanzi et al., 2016;Magnusson et al., 2015).This explains why hydrology-oriented cryospheric models tend to have low complexity when it comes to internal layering and micro-scale properties.Fourth, processes relevant to cryosphere water resources span horizons from a few hours (such as rain-on-snow events, see Würzer et al., 2017) to decades (such as glacier dynamics, see Huss et al., 2010), implying that models used for real-world forecasting must be efficient enough to provide landscape-scale predictions in a timely manner (say, a few hours, see Pagano et al., 2014).Ultimately, these four factors trace back to empiricism rather than reductionism being the dominant (and perhaps most successful) paradigm in hydrology (Savenije, 2009), owing to unresolved issues related to upscaling mechanicistic laws to the landscape and measuring the complete heterogeneity of hydrologic processes (Blöschl and Sivapalan, 1995;Beven, 2006).
Here, we present Snow Multidata Mapping and Modeling (S3M) v5.1, a snow and glacier model developed by CIMA Research Foundation (https://www.cimafoundation.org/).Specific aspects of interest in S3M are (1) a spatially explicit prediction of both dry and wet-snow spatial patterns, as well as bulk snowpack liquid water content (θ W , in vol%), an increasingly decisive variable for snowmelt and avalanche hazard forecasting (Techel and Pielmeier, 2011;Mitterer et al., 2013;Wever et al., 2014;Avanzi et al., 2015;Wever et al., 2016); (2) the combination of both snow and glacier mass dynamics in a coherent modeling framework, including the so-called ∆h parametrization by Huss et al. (2010) and melt beneath supraglacial debris; (3) provisions for assimilating various decision-relevant variables like SWE, snow depth, and satellite-based snow-cover area.S3M v5.1 is the last generation of a model originally proposed by Boni et al. (2010), but significantly developed thereafter (henceforth, simply S3M).
The paper is organized as follows: Section 2 focuses on model description, including both snow and glaciers.Section 3 presents an example of results for an inner Alpine valley in north-western Italy where various versions of this model have been operational since the 2000s (Aosta valley).Finally, sections 4 discusses model applicability and future developments.The Supporting Information includes an User Manual discussing run preparation, execution, and post-processing.
The main modeling philosophy behind S3M is to provide a ready-to-use tool for applications over large areas and when a short turnaround is needed.Thus, S3M complies with the four pragmatic factors of operational, hydrology-oriented cryospheric models outlined in the Introduction, including being spatially distributed, parsimonious as for both input-data requirements and complexity, and enough computationally efficient to be deployed in operational, real-time flood-forecasting chains (Laiolo et al., 2014).These chains generally ingest an ensemble of weather predictions, may include parameter and/or state perturbation, assimilate remote-sensed and in-situ measurements, and must provide predictions for a variety of closure sections across the landscape in a matter of hours, if not minutes (Pagano et al., 2014).As such, they tend to be much more conceptual and computationally simple than research-oriented Earth-system models (Pagano et al., 2014).S3M tries to bridge the gap between these two realms, by accompanying simplicity and computational efficiency with an open-source environment that allows for frequent updates and improvements.
Parsimony implies a number of trade-off choices regarding what process representation to include, and to which extra cost.
Current model physics include precipitation-phase partitioning, snow and glacier mass balances, a hybrid temperature-index and radiation-driven melt parametrization, snow rheology and hydraulics, and a data-assimilation protocol.Relevant drivers of snowpack and glacier evolution that are not yet included comprise internal snowpack energetics and the full energy balance, longwave losses, turbulent heat fluxes, sublimation, and canopy-snow interactions.While recent literature has showed that the added value of complex, physics-based snow models over more parsimonious alternatives for variables that are relevant to hydrology may sometimes be elusive (Rutter et al., 2009;Magnusson et al., 2015;Zaramella et al., 2019;Girons Lopez et al., 2020;Günther et al., 2019), designing more holistic and physics-based operational models remains a key to achieving the most accurate representation of snow temporal dynamics and spatial patterns, especially at fine resolutions (see for example results in Lafaysse et al., 2017;Vionnet et al., 2021).Strategies to include these processes in future releseases of S3M are discussed in Section 4.
S3M is a raster-based model, with the same set of equations being solved for each cell and no lateral mass or energy transfer besides glacier change in thickness.All input, state, and output variables are distributed, meaning they are passed to the model as rasters with fixed resolution in geographic degrees (see the Supporting Information).Spatial resolution can be set by the user (we will consider a ∼240 m resolution in our example of application in Section 3).All equations are ordinary differential equations with no need for iterative computations, so they are solved for all pixels in the simulation domain using a forward-Euler method, as this method provides comparatively high numerical stability and minimizes computational time.The time step of the model is flexible, but it is generally set to 1 hour.
We define snow and glacier ice as a mixture of three constituents: ice, liquid water, and air.Following De Michele et al. (2013) and Avanzi et al. (2015), the control volume of unitary area (h T OT ) for each pixel of the simulation domain is defined as: where h G is glacier thickness (m), h S is the height of snow (often referred to as snow depth, m), V tot is the control volume in m 3 and A is the area of the pixel (m 2 ); h G = V G /A and h S = V S /A, where V G and V S are the total volume of glacier and snow within the pixel under study.We define M S and M G as the mass of seasonal snow and glacier ice for each pixel, respectively (in kg).As for seasonal snow, this mass is Snow is a foam of ice (Kirchner et al., 2001), meaning that the dry constituent occupies a porous skeleton of height h D (m, volume V D ) and porosity n = h P /h D (-), where h P is the height of pores in the mixture (m, volume V P ).Thus, we define the density of the solid skeleton ρ D (in kg m −3 ) as: where ρ i = 917 kg m −3 is ice density.The density of the wet constituent is equal to that of liquid water: 1000 kg m −3 (V W is the volume of liquid water in the mixture, height h W ). The density of glacier ice is assumed equal to ρ i , with no progressive compression and/or air expulsion from cavities (Paterson, 1994).
During most of the snow season, the dry and wet constituents of the snowpack are both contained in h D , the prevalent volume.However, an oversaturation condition takes place during the last instants of the snow season due to phase change (De Michele et al., 2013).Despite being a limit, virtually unmeasureable scenario, including this oversaturation condition is important from a numerical-stability standpoint.Thus, the total control volume of snow h S (that is, snow depth) is hereby defined as: where are Macaulay brackets, which provide the argument if this is positive, otherwise 0. In other words, h S is equal to the height of the porous structure h D plus -if present -the oversaturated volume h W − nh D (De Michele et al., 2013).
Accordingly, the bulk snow density ρ S (in kg m −3 ) is and Snow Water Equivalent (in m w.e.) is We also define the bulk volumetric liquid water content of the snowpack as: Glaciers are modeled as a single-phase material, the mass of liquid water and air being being negligible compared to glacier ice.Thus, the Ice Water Equivalent for glaciers (IW E) is: Figure 1 summarizes the prognostic variables and main mass fluxes of S3M.Dynamic inputs for S3M are total precipitation, air temperature, relative humidity, and incoming shortwave radiation (see the Supporting Information for details).

Snow: mass-conservation equations
The mass-conservation equations for the dry and wet constituents of the snowpack read as follows: where Ŝf is the snowfall mass flux, M is the snowmelt mass flux, R is the refreezing mass flux, Rf is the rainfall mass flux, and Ô is the outflow mass flux (also known as snowpack runoff, see Avanzi et al., 2019).All these mass fluxes are expressed in kg ∆ t −1 and are denoted with aˆ, as opposed to customary hydrologic fluxes in mm ∆ t −1 , which will be denoted withoutî n the following.Given that M D = ρ D h D A, we simplify Equation 7as hG is glacier thickness (m), MG is ice melt (mm ∆t −1 ).Note that albedo is here listed as prognostic variable for clarity, but the actual, underlying prognostic variable is technically snow age As.The background image is Rutor glacier in north-western Italy.
where SW E D is the dry-snow water equivalent (ρ D h D ρ −1 W ) and S f , M , and R are the snowfall, snowmelt, and refreezing mass fluxes in mm w.e.∆t −1 .Note that SW E D and all related mass fluxes will henceforth be expressed in mm w.e., with a conversion by 1000 mm m −1 being implicitly included between Equation 7 and 10.Likewise, we simplify Equation 8 as with SW E W being the wet-snow water equivalent, and R f and O being the rainfall and snowpack-runoff mass flux in mm w.e (again, note that SW E W and all related mass fluxes will henceforth be expressed in mm w.e., with a conversion by 1000 mm m −1 being implicitly included between Equation 8 and 12).
Equations 10 and 12 are the two fundamental snow mass-conservation equations of S3M, which thus offers a spatially explicit, prognostic simulation of both the dry and wet constituents of snow (total SW E being instead a diagnostic variable: S3M assumes that simulated SWE, as well as all other model variables and fluxes, are representative of spatially averaged snowpack conditions across the simulated pixel.Whilst this assumption is coherent with current practice in many hydrologic models, it does not take advantage of scaling mechanisms with fractional snow-covered area to fully reconstruct sub-grid snow-depletion curves.

Snow: mass-flux parametrizations
Mass fluxes in Equations 10 and 12 requiring specific parametrizations are snowfall and rainfall (R f and S f ), snowmelt and refreezing (M and R), and snowpack runoff (O).

Precipitation-phase partitioning
Snowfall and rainfall in S3M are estimated from total precipitation (P ), an input for the model; precipitation-phase partitioning is based on the empirical approach described in Froidurot et al. (2014): where p s and p r are the probabilities of snowfall and rainfall, respectively, α, β, and γ are fixed parameters derived by Froidurot et al. (2014), T air is air temperature in °C, and RH is relative humidity in %.S3M assumes p s and p r to be equal to the actual proportions of rainfall and snowfall over total precipitation.Following Froidurot et al. (2014) and references therein, α = 22, β = −2.7,and γ = −0.2.

Snowmelt and refreezing
Snowmelt (M ) is computed if both concurrent T air and mean air-temperature over the previous 10 days ( T10d ) are greater than, or equal to, T τ , a user-defined threshold usually assumed equal to 1 °C (Pellicciotti et al., 2005); otherwise, M = 0.
The first condition is standard in degree-day models and accounts for snow melt occurring during periods with a supposedly positive energy balance (meaning a net gain of energy for the snowpack).The second condition is a novel addition of S3M to the literature to keep track of cold content using a parsimonious approach (cold content being a measure of the snowpackenergy deficit to be satisfied for actual melt to start, see Jennings et al., 2018).The basic idea is that T10d evolves with a certain delay compared to T air , so that setting an additional threshold on T10d helps avoiding non-physical melt during short warm spells that come after a somewhat long cold period.Other simple approaches to estimate cold content exist, such as that based on mean-seasonal temperature by Schaefli and Huss (2011), but are also in their early stages.Our approach is the result of intensive trial and error and specifically aims at suppressing erroneous mid-winter melt episodes that do not appear in validation data (see Section 3).
Snowmelt is parametrized using a hybrid physics-based and degree-day approach decoupling radiative forcing from temperaturedriven melt (similarly to Pellicciotti et al., 2005): where S r is incoming shortwave radiation (an input for S3M, in W m −2 ), α S is snow broadband albedo (-), λ f is the specific latent heat of fusion (0.334 MJ kg −1 ), m r is a degree-day parameter (mm °C−1 day −1 ), m rad is an adimensional modulating factor to convert shortwave radiation into actual melt (similar to an efficiency parameter, see below), and c M is a timestepadjusting parameter.While S r is internally converted to MJ m −2 according to the model timestep (see ∆t in Equation 15), the air-temperature part of Equation 15 is computed by first considering an equivalent day with average temperature equal to T air , regardless of the model time step.The snowmelt part depending on air temperature is then readjusted by c M = ∆t/86400 to pass from one day to the actual time step (∆t is in seconds).This workaround is due to the degree-day approach being originally conceived for daily applications (see Hock, 1999Hock, , 2003;;De Michele et al., 2013, and references therein).Note that S3M internally sets S r to 0 between 7PM and 7AM, in order to remove spurious noise in radiation sensors (an alternative approach being the computation of sunrise and sunset time).Proper unit conversions are implicitly included in the radiation part of Equation 15to first pass from J m −2 to MJ m −2 and then from m w.e. to mm w.e.
Broadband albedo is computed once per day at midnight according to Laramie and Schaake (1972): where τ α is an albedo-decay coefficient equal to 0.12 d −1 if average air temperature over the previous 24 hours is higher than 0°C, and otherwise equal to 0.05 d −1 .A s is snow age (in d), defined as the number of days since the last significant snowfall.S3M considers as significant snowfall one day with at least 3 mm of total snowfall.Snow age A s is updated every day at midnight: if cumulative snowfall during the previous 24 hours is less than 3 mm, then snow age is increased by 1 day; if not, than snow age is reset to 0.
Similarly to other snowmelt models (Rango and Martinec, 1995), melt parameters m r (mm °C−1 d −1 ) and m rad (-) are calibration-based, although the sensitivity of S3M to both is rather low (see results in Section 3).This low sensitivity is likely because explicitly separating the radiation-and temperature-driven components of melt brings these parameters closer to a fully first-principals energy balance model than standard degree-day approaches.
An explicit separation between the temperature-and the radiation-driven components has been underused in the literature.
For example, the seminal work by Hock (1999) did include potential radiation, but embedded it into a degree-day parameter rather than explicitly separating the radiation and the temperature components.Follum et al. (2015) replaced the temperature term with a proxy from a radiation balance, which may be suitable in regions where snowmelt is mostly radiation driven (e.g., the western US, see Bales et al., 2006), but would need some form of temperature dependency in temperate regions like the Alps.Pellicciotti et al. (2005) are among the few examples where the radiation and the temperature component are fully decoupled, but they focused on glacier ice during summer.This isothermal condition is very efficient for shortwave radiation to convert into actual melt.However, applying the original approach by Pellicciotti et al. (2005) to snow revealed a tendency to overestimate melt rate early in the snowmelt season, because it assumes net shortwave radiation to translate into melt regardless of the actual cold content.In subfreezing conditions, in fact, a fraction of net shortwave radiation is used to raise snow temperature, a mechanism that becomes increasingly unimportant as the season progresses and snow conditions tend towards isothermal.
To mimic this transition from subfreezing to isothermal conditions, we propose a modification to the original approach by Pellicciotti et al. (2005) in the form of a novel 10-day-temperature-modulated efficiency parameter m rad that increases with T10d according to a sigmoid function (Figure 2, a): While parameter m rad is user-defined, we note that m rad ∼ 1.10 means that m rad → 1 when T10d > 10 °C.This corresponds to a 1:1 conversion of net shortwave radiation into melt when isothermal conditions like those by Pellicciotti et al. (2005) dominate.On the other hand, m rad → 0 when T10d → 0 °C.m rad is set to 0 if the equation above predicts a negative value.Note that this temperature-modulated efficiency parameter is a proxy of cold content, but it does not imply an explicit computation of cold content.It is a first attempt to take into account thermal inertia in subfreezing conditions and how it is related to external climate, similarly to Schaefli and Huss (2011).At the present stage, no relation with snow depth or internal snow temperature is included.
S3M considers a similar relation between T10d and m r through a sigmoid function and a user-defined tuning parameter (m r , see Figure 2, b): Here again, we note that m r ∼ 1.40 mm C −1 d −1 corresponds to ∼0.05 mm C −1 h −1 when T10d > 10 °C, which agrees with estimates by Pellicciotti et al. (2005) in isothermal conditions on ice.While establishing a relationship between melt parameters and T10d is novel, previous literature has already suggested a seasonal variability in the degree-day parameter that is conceptually similar to our approach (see Bongio et al., 2016, and references therein).Also, note that m r should be much smaller than degree-day parameters listed by Hock (2003), because the latter were supposed to account for both the temperature-driven and radiation-driven component of snowmelt.Parameter m r is set to 0 in case the equation above returned a value lower than 0.
Refreezing (R) is computed when T air < T τ , using a simple degree-day approach as in Avanzi et al. (2015): Differently from Avanzi et al. (2015) or Schaefli et al. (2014), we do not decrease m r by a reduction factor when computing refreezing.In standard degree-day models that do not separate the temperature and radiation components of melt, that reduction factor accounts for refreezing melt rate being smaller than snowmelt rate for the same temperature difference, given the usual lack of incoming-shortwave radiation during refreezing-prone periods.This reduction factor is not necessary in S3M, because the contribution of incoming-shortwave radiation is explicitly accounted for in Equation 15 and excluded in Equation 19; in other words, m r only accounts for turbulent and longwave-radiation factors.

Snowpack runoff
While the term runoff in catchment hydrology generally denotes overland flow, we follow here customary nomenclature in snow hydrology and call snowpack runoff the amount of liquid water discharged by the snowpack (Wever et al., 2014).This flux is parametrized according to a matrix-flow approximation (Colbeck, 1972;Avanzi et al., 2015Avanzi et al., , 2016)): where µ W is water dynamic viscosity, g is acceleration due to gravity, K W is the intrinsic permeability of water in snow (m 2 ), and α is a time-and unit-conversion parameter (equal to 1000 mm m −1 ×∆t, with ∆t in seconds).Assuming ρ W = 1000 280 kg m −3 , g = 9.81 m s −2 , and µ W = 1.79 × 10 −2 kg m −1 s −1 (DeWalle and Rango, 2011), then with α = 5.47 × 10 5 m −1 s −1 .We predict K W following again Colbeck (1972): where K is the intrinsic permeability of snow (in m 2 ) and S is the effective saturation degree: Sr i is the irreducible saturation degree computed based on Kelleners et al. (2009) D is the saturation degree.Snowpack runoff is set to 0 if Sr < Sr i .Intrinsic permeability of snow is predicted based on Calonne et al. (2012): where r e is the equivalent sphere radius (m), a conceptual, characteristic length of snow microstructure corresponding "to the radius of a monodisperse collection of spheres having the same specific surface area (SSA) as the sample considered" (Calonne et al., 2012).Variable r e or SSA (m 2 kg −1 ) are likely the most objective metrics of snow microstructure to date (Carmagnola et al., 2014), traditional grain size being subjective and cumbersome to measure (Fierz et al., 2009).Still, SSA and r e are insufficient to fully characterize snow structure, because grain shape also plays an important role in the two-point correlation function (Krol and Löwe, 2016).We relate r e to SSA by definition: and diagnostically estimate SSA following Domine et al. (2007): with SSA being SSA in cm 2 g −1 and ρ D is ρ D in g cm −3 .We used Equation 26 to predict SSA because S3M does not include a prognostic simulation of snow microstructure.However, Domine et al. (2007) clearly show that density alone is a modest predictor of SSA (R 2 = 0.43).
In order to avoid high saturation values in a shallow snowpack to cause large outflow rates and non-physical negative values of SW E W , Equation 20 is used as long as Sr < 0.5 or SW E D > 10 mm w.e..In situations with Sr ≥ 0.5 or SW E D < 10 mm, Equation 20 is bypassed and SW E W is directly converted into snowpack runoff.

Dry-density equation
Differently from liquid water, bulk dry-snow density in S3M is not invariant with time because of three main factors: compaction, new-snow events, and refreezing: These factors do not include snow metamorphisms (Pinzer, 2009), mainly because of the scale mismatch between such processes and the one-layer approach of S3M.
Regarding compaction, we assume a linear profile for stress with depth and start from the momentum-conservation equation for a representative element at 66% depth, which experiences an average stress: with σ v being vertical stress.Equation 28 is then coupled with a viscous rheological equation to obtain, via the definition of the vertical strain rate: with η being viscosity.We finally follow De Michele et al. ( 2013) and references therein and define viscosity as an exponential function of dry-snow density and snow temperature: or in a more compact form as: T S is snow mean temperature (°C), c 1 = 0.001 m 2 h −1 kg −1 , and c ∆t is a timestep-adjusting coefficient (∆t × 3600 −1 s −1 h, with ∆t in seconds).
Because S3M does not solve for the full energy balance, it also does not simulate snow-temperature profiles.The fact that snow temperature in S3M has no implication for snow melt, compounded by the sensitivity of the settling equation to snow temperature being rather small (not reported for brevity) led us to introduce a simple parametrization here: if T air < 0°C, snow mean temperature is assumed to follow a linear profile between snow-surface temperature and ground surface (assumed equal to T air and 0 °C, respectively, following pieces of experimental evidence like Ohara and Kavvas, 2006;Filippa et al., 2014), while otherwise it is set to 0°C.
New events change bulk-snow density proportionally to snowfall depth versus existing-snow depth: We handle refreezing with a similar approach to new events:

Data assimilation
The assimilation framework of S3M is a result of CIMA's operational-forecasting procedures as summarized into the Flood-PROOFS suite (Laiolo et al., 2014;Avanzi et al., 2021).These procedures -external to S3M -include generating maps of snow depth and satellite-based scene classification (hence, snow-covered area -SCA), as well as processing SWE maps from third parties (e.g., from interpolation of ground manual measurements, see Avanzi et al., 2021).Given that snow depth and satellite maps are generated with a comparatively high frequency (up to daily), their assimilation in S3M is performed in correspondence to the timestamp to which they refer (this nominal timestamps must be the same for both snow-depth and satellite maps, collectively referred to as Updating maps).SWE maps have various temporal frequencies (usually, weekly), thus S3M allows the user to specify a temporal window of influence, that is, a period after the official issue date of the SWE map during which the map is assumed to be valid.
Assimilation of Updating maps additionally requires a Kernel map and a Quality map.The Kernel (K) is generally an output of the geostatistics-based interpolation method employed to generate the snow-depth maps and is used to optimally combine observations and model predictions (see below).Instead, the Quality map is used to automatically skip assimilation when values in this map are below an user-defined threshold.For example, the operational convention in Flood-PROOFS is to forego assimilation during days with large cloud obstruction.Thus, in Flood-PROOFS the Quality map for a given day is computed as the ratio between pixels classified as snow or ground over the total number of pixels (so, in fact, this map reports the same scalar for each pixel).In this way, quality is a measure of the proportion of the satellite map that is covered by clouds.We stress, however, that this quality can be defined based on any user's need, with S3M skipping assimilation below a quality threshold even if the corresponding snow-depth map was available.
Prior to assimilation, S3M blends information from the snow depth and the SCA maps based on quality.For example, K is doubled wherever quality is 2.5 times the quality threshold and the satellite indicates ground (ID = 0).Also, snow-depth maps are set to 0 wherever quality is greater than the quality threshold and the satellite indicates ground.Finally, snow-depth maps are set to missing values wherever quality is lower than the quality threshold, modeled SWE is less than 20 mm w.e., and the satellite indicates clouds, no decision, or no data (ID = 2, 3, and -1, respectively).
SW E maps, on the other hand, use neither a quality flag nor a spatially distributed Kernel, again owing to how these maps are received and handled by Flood-PROOFS.For assimilating SWE maps, K is thus assumed to be a function of time elapsed since the issue date, with no spatial variability: where t SW E is the official issue time of the SWE map, t is time relative to this issue date (in days), and σ is equal to half of the validity days after SWE-map issue date (user-defined parameter, see the Supporting Information).On the same date of the issue date, t = t SW E and so K = W , where W is a user-assigned maximum weight of the map to be assimilated.No quality threshold is used for SWE maps given that they are usually the result of reanalysis rather than real-time automatic processing.
Currently, S3M performs data-assimilation exclusively in the form of SWE.Thus, snow-depth-map assimilation requires a preliminary step to convert snow depth into SWE via modeled bulk-snow density: where h S,obs and ρ S,S3M are observed snow depth according to the snow-depth map and simulated bulk-snow density according to S3M, respectively.This step implicitly assumes that snow density is a less relevant source of uncertainty than snow depth in estimating SWE, which is supported by snow-density temporal patterns being consistent from year to year (Mizukami and Perica, 2008).
Updating and SWE maps are assimilated into S3M using a Newtonian-relaxation approach: where SW E S3M,post and SW E S3M,prior are the a-posteriori and a-priori SWE.Note that Newtonian Relaxation (also known as nudging) is different from direct insertion, where model estimates are directly replaced by observations; in a nudging scheme, the correction factor is proportional to the difference between observations and model outputs via a Kernel weight (Boni et al., 2010;Mazzoleni et al., 2018).
After assimilating bulk SWE, a few prognostic variables of S3M are modified through factor U SW E : This step is needed since total SWE in S3M v5.1 is only a diagnostic variable, and assimilating it does not affect model predictions unless the true prognostic variables are also modified (in this case, SW E D and SW E W ). Factor U SW E assumes that both the dry and the wet phases are proportionally affected by data assimilation.It is also assumed that dry-snow density does not change during assimilation.Given that dry-density evolution does depend on SW E D , this is a simplification.
S3M optionally supports assimilating only positive differences in Equation 36, that is, only correcting modeled SWE if observations are larger than simulations.This experimental configuration helps when assimilating observed SWE aims at correcting for precipitation under-catch (Ryan et al., 2008).Such an approach is, e.g., standard in avalanche-forecasting models like SNOWPACK (Lehning et al., 2002).However, assimilating only positive differences will override the SCA component of the assimilation package, because SCA assimilation in S3M is performed indirectly via setting observed snow depth to zero in areas with observed ground.
The last step in the snow component is to perform a set of sanity checks (e.g., set to zero all state variables where SWE→ 0) and to compute the snow mask, a binary map with the same size of the simulation domain reporting 1 where SW E > 0.1 mm, and 0 elsewhere.This mask is an important output of S3M that is sometimes used by hydrologic models to adjust process representation in areas of snow (for example, inhibiting evapotranspiration).

Glacier component
The glacier component of S3M offers three alternative modules: (1) a simple, melt-only approach with no mass balance and no snow-to-ice conversion (G1), which is usually the default choice for short-term flood-forecasting-oriented simulations; (2) a melt-only approach with mass balance but no parametrization of glacier dynamics and no snow-to-ice conversion (G2); (3) an approach with a full mass balance, a parametrization of glacier dynamics (the ∆h parametrization), and snow-to-ice conversion (G3).The last two approaches are the most suited options for long-term, climate-scenario-oriented simulations.

G1: Melt-only approach
In the most basic approach (G1), glacier melt takes place on snow-free glacier pixels with a similar parametrization as that of snow (Equation 15), with only two changes.The first is that glacier albedo α G is constant (0.4, see Davaze et al., 2018), as opposed to snow albedo changing with snow age.The second change is that ice-melt rate (M G ) on debris-covered glaciers may be corrected compared to the theoretical melt rate on debris-free glaciers (M G ) using a multiplicative reduction factor, f debris (Huss and Fischer, 2016): The correction factor f debris can be estimated with various approaches, for example following Huss and Fischer (2016) who prescribed it based on the so-called Østrem curve.Accordingly, values of f debris are generally smaller than 0 for very shallow debris (up to ∼5 cm), and between 0 and 1 otherwise (Nicholson and Benn, 2006).S3M expects f debris to be a spatially distributed, input parameter included in the so-called static-data suite (see the Supporting Information).Glacier melt is directly converted into ice runoff, with no routing (Figure 1).Glacier pixels are defined based on a glacier mask (see the Supporting Information).

G2: Melt-only approach with mass balance
This approach (G2) is equivalent to G1 with the only difference that glacier thickness (h G ) is dynamically updated every hour based on glacier melt.To do so, ice melt in mm w.e. according to Equation 15 and optionally Equation 40 is converted to meters of ice using ice density ρ i .The mass-conservation equation reads: Wherever h G → 0 as a result of multi-year melt, ice melt on that pixel is not computed anymore.S3M expects h G to be a spatially distributed input (included either in the so-called restart data or in the static data, see the Supporting Information).
Spatially explicit datasets of h G could come from either in-situ surveys (e.g., see Colombero et al., 2019) or from estimates based on the surface mass balance and ice-flow velocities (e.g., see Rabatel et al., 2018).

G3: mass-balance approach with glacier dynamics and snow-to-ice conversion
This approach -G3, ideal for multi-year simulations -includes snow-to-ice conversion and a specific mass-redistribution approach called the ∆h parametrization (Huss et al., 2010), which allows one to implicitly account for glacier-movement effects without implementing a full ice-flow model.Because the ∆h parametrization is better suited for multi-year time scales rather than day-to-day thickness changes, module G3 requires the user to define the month of water-year start, so that S3M will accumulate glacier melt for each pixel throughout the water year and update h G at the beginning of each new water year.
Regardless of this accumulation procedure, ice melt is still outputed every time step, so that seasonality in runoff-generation processes is preserved.In other words, this accumulation procedure only regards changes in glacier thickness and not ice-runoff generation.This is important in case S3M was coupled with a hydrologic model.
Snow-to-ice conversion is performed by simply prescribing that, on pixels with h G > 0, any residual SWE at the end of each water year is added to h G .Consequently, SWE as well as all snow-related bulk state variables are reset to 0. This reset is not performed in areas where SWE> 0 and h G = 0 at the end of the season; in such conditions, the snowpack is maintained through the start of the new water year.This approach is less rigorous than considering firn, the intermediate step between snow and ice, even though previous work by Schaefli et al. (2005) in Switzerland showed that the use of a separate degree day factor for firn may not significantly improve hydrologic predictions.
As for the ∆h parametrization, this is presented in Huss et al. (2010) and further discussed for hydrologic models by Seibert et al. (2018), so we limit ourselves to a short overview here.This approach starts from the empirical intuition that glacier-thickness changes as a result of both the mass balance and glacier flow have recurring patterns throughout seasons of persistently negative mass balances (Huss et al., 2010).By parameterizing these recurring patterns, the ∆h parametrization allows one to simulate the effect of movement in addition to the mass balance without a complex ice-flow model.These patterns are derived through differentiating two digital surface models (Huss et al., 2010) and then fitting a glacier-specific power law like where a, γ, b, c are calibration parameters, ∆h is the change in surface elevation (normalized by the maximum decrease across all glacier pixels), and h r is normalized glacier elevation defined as with h max and h min being the maximum and minimum elevations of that glacier at the beginning of the water year and h being glacier elevation at a given pixel, respectively.Note here that h, h max , and h min are elevations, not thicknesses like h G .
In practice, applying the ∆h parametrization requires (1) assigning an ID to each glacier for which the user would like to use the ∆h parametrization (S3M expects this ID to be a positive integer); (2) mapping these glaciers ID on the simulation raster, so that S3M will be able to identify all pixels of the simulation domain belonging to a given glacier; (3) a-priori deriving Equation 42for each glacier of interest and passing it to S3M as a pivot table, where ∆h is sampled for a number of discretized h r .S3M will then assign a ∆h for each pixel of a given glacier using a nearest-neighbor interpolation of this pivot table.Both the glacier-ID map and the pivot table are part of the so-called static-data suite in input to S3M (Supporting Information).
Once these preliminary steps are performed, S3M computes ice melt as in module G2, but h G is not dynamically updated.
Instead, ice melt for each pixel is accumulated to yield b a , the cumulative mass balance in mm w.e.At the end of the water year, this information is used to compute factor f s , which is employed to scale Equation 42 (the ∆h profile) and so derive the actual change in glacier thickness for each glacier pixel: where i denotes the i-th pixel of a given glacier, N G is the total number of pixels of that glacier, A i is the area of each pixel, and A tot is the total area of the glacier.Once h G,i → 0, that pixel is excluded from all computations and h max and h min as well as all other variables are updated accordingly.
The general mass-conservation equation for glacier pixel i using G3 therefore reads: with SW E i,r the residual snow water equivalent on that pixel at the end of the water year, while IW E i,W Y +1 and IW E i,W Y are ice-water equivalent for pixel i during the current and the previous water year (W Y ), respectively.
The implementation of the ∆h parametrization in S3M provides a number of modeling degrees of freedom.First, S3M assumes a non-dynamic mass balance for all pixels that are part of the model glacier mask, but have no glacier ID assigned; the user can explicitly choose this approach also for glaciers having a glacier ID, by setting all entries of the corresponding pivot table to -9999.This option is useful either where fitting a ∆h parametrization is cumbersome (e.g., spatially incoherent glaciers) or for glacial remnants that are not moving anymore (e.g., glacierets).Second, the modeler can use different ∆h parametrizations for various parts of the same large glacier, by simply assigning different glacier IDs to these parts and providing specific pivot-table entries for each ID.Also note that the ∆h parametrization is bypassed whenever one glacier occupies only one pixel.
3 Case study: Aosta valley, NW Italian Alps S3M is open software, including algorithms to prepare input data and set up computational environments and libraries.Links to all code are reported in the User Manual (see the Supporting Information), with a general reference being CIMA Foundation's Hydrology and Hydraulics repository at https://github.com/c-hydro(Github organization).
This Section presents an application of S3M for an inner-Alpine valley located in north-western Italy (Aosta valley, Figure 3).This area has steep elevation gradients, with the main valley at elevations on the order of 300-400 m ASL and peaks as high as 4800 m ASL (Mont Blanc) or 4478 m ASL (Matterhorn).About 4% of Aosta valley is covered by glaciers (134 km 2 ), some of which are characterized by thick debris and extend over several kilometers (e.g., the Miage glacier in the Mont-Blanc massif).With its cryosphere-dominated water supply and complex precipitation-topograghy interactions leading to marked rain shadows (Avanzi et al., 2021), Aosta valley is a formidable test bed for S3M.
S3M has been operational in Aosta valley since the early 2000s, as a component of a flood-forecasting chain called Flood-PROOFS (see Laiolo et al., 2014;Avanzi et al., 2021, and Section 2).This chain includes algorithms to spatialize and downscale weather-input data of precipitation, air temperature, relative humidity, and radiation, automatically generate daily maps of snow depth and use MODIS snow-covered area, and process independently derived weekly maps of SWE (see the Supporting Information).Together with runs of S3M in assimilation mode at ∼240 m spatial resolution, these tools inform real-time forecasts of streamflow at relevant closure sections (Laiolo et al., 2014;Avanzi et al., 2021).Details about spatialization techniques and hydrologic modeling in Flood-PROOFS are available elsewhere (Boni et al., 2010;Laiolo et al., 2014;Avanzi et al., 2021) and are not discussed here for brevity.In the present paper, we instead leverage our application in Aosta valley to provide guidelines on how to calibrate S3M in a real-world case study (Section 3.1) and how to validate and interpret model results for the snow (Sections 3.2 to 3.4) and the glacier component (Section 3.5).

Calibrating S3M
S3M parameters that can be calibrated include the radiation and temperature snowmelt parameters m rad and m r , the threshold temperature for inhibiting snowmelt (T τ ), the maximum weight of SWE-assimilation maps (W ), as well as a number of climatological thresholds used to constrain model predictions (e.g., maximum and minimum snow density, or the snow-quality threshold to enable data assimilation, see the Supporting Information).While previous work from the author team in Aosta valley has involved calibration of many of the above parameters, the two melt parameters, m rad and m r , remain the prime 505 calibration parameters of this model -like many other snowmelt models (Hock, 2003;Pellicciotti et al., 2005).
Because S3M is employed in assimilation mode in Aosta valley, and this mode includes MODIS snow-covered area, our calibration rationale focused on maximizing fit for point predictions rather than for spatial patterns.Still, calibration was performed considering open-loop simulations to avoid model performance to be spuriously driven by assimilation rather than parameter values.Calibration data comprised spatially distributed continuous-time measurements of 53 snow-depth sensors and temporally discontinuous manual measurements of snow depth collected across the region for water-supply or avalanche forecasting (see Avanzi et al., 2021, and Figure 3 for a data inventory).The calibration period covered water years 2010 through 2019, where the bulk of data was concentrated; the water year is a period between September 1 to August 31 and is indicated with the calendar year when it ends.
Our calibration protocol was based on performing multi-year simulations of S3M for a range of values of m rad and m r : [0.8, 2] and [0.5, 1.5] mm °C−1 day −1 , respectively.These ranges were chosen based on the fact that m rad ∼ 1.1 implies a 100% efficiency of transmitted shortwave radiation in generating melt under likely isothermal conditions ( T10d → 10 • C, see Figure 2(a) and Section 2), while m r ∼ 1.4 tallies with previous work by Pellicciotti et al. (2005).The parameter space was explored for increments of 0.025 of m rad and 0.05 mm °C−1 day −1 of m r , and we first calibrated an optimal value of m rad and then of m r .This meant running 29 decade-long simulations for all tentative values of m rad , finding the optimal one, and then re-running 25 simulations tuning m r .This sequential calibration approach is only one out of several possible calibration approaches, including strategically exploring the parameter space (Razavi et al., 2019).Here, we chose a sequential approach for this illustrative case study mainly for computational-resource constraints.
As objective metric, we minimized Obj(Θ) = 0.5 × (1 − KGE(Θ)) + 0.5 × (RMSE(Θ))/RMSE(Θ), with Θ being the generic parameter value, KGE being the Kling-Gupta Efficiency metric as defined in Kling et al. (2012), RMSE being the Root Mean Square Error, and RMSE(Θ) being the mean RMSE across all parameter values.This objective metric was chosen to combine the features of KGE (and in particular its focus on correlation, bias, and ratio of coefficients of variation) with a specific weight for large errors (RMSE).The RMSE was normalized by its mean across all parameter values to make it adimensional and so comparable to KGE.For each parameter value Θ, Obj was computed between observed snow depth (Figure 3) and simulated h S , both merging all data from all years in one sample and separately for each water year (all-years and yearly values, respectively).
Median yearly values across all water years show a minimum for m rad = 1.125 and then m r = 1.10 mm °C−1 day −1 , in line with expectations (Figure 4 and Section 2).Objective-metric values showed remarkable variability across water years (see the quartile range in Figure 4), owing to significant variability in the original sample of snow-depth values between warm and cold years (Avanzi et al., 2021).Thus, multi-year calibration periods are recommended, although the range of variability in median and all-year Obj in Figure 4 suggests that the sensitivity of S3M to m rad and m r is low, especially if one considers that calibration-based snow models are prone to large drops in performance outside the calibration sample (Hock, 2003;Avanzi et al., 2016).We interpret this to be because of the hybrid physics-based and temperature-index approach employed in S3M to predict snowmelt (Equation 15), which appears to better constrain parameter values than a purely temperature-index approach because it brings S3M closer to a first-principals energy balance model.

Evaluation: point snow depth
Figure 5 show simulated vs. observed snow depth (symbol HS per guidelines by Fierz et al., 2009) for the 52 snow-depth sensors in Figure 3, the snow-depth sensor at the Torgnon study plot, as well as for all manual snow-depth measurements taken for avalanche, hydropower, and water-supply forecasting (panels a to e, respectively).The evaluation period in Figure 5 is water years from 2004 to 2009 and 2020, that is, all years before and after the calibration period in Section 3.1.Simulations in this and following sections were carried out in assimilation mode, as this is the approach we generally use in real-time forecasting and so results are more representative of real-world model performance.More details on this assimilation procedure are reported in Avanzi et al. (2021); note that assimilation is not performed between May and July to avoid interferences with the simulation of the depletion phase of the seasonal snowpack.
Snow-depth-sensor and weekly-snow-sample measurements in Figure 5(a) and (e) were indirectly assimilated in S3M as they are involved in the computation of the so-called Updating and weekly-SWE maps (see Section 2), meaning their performance statistics are only reported for reference.Non-assimilated data at Torgnon (Figure 5(b)) and from avalanche probes (Figure 5(c)) maintain comparatively high values of RMSE and KGE (12-37 cm and 0.66-0.72,respectively), with no evident tendency for systematic under-or overestimations.Open-loop results were not reported here for brevity, but generally show comparable performance for these datasets.Thus, we conclude that our calibration strategy not only showed little sensitivity of the model to m rad and m r (Figure 4), but also led to a robust performance across nearly twenty water years and for areas that were not included in the calibration pool.
Peak-snow-depth courses show a significantly lower performance, in particular in terms of a larger dispersion around the 1:1 line, a larger RMSE (122 cm), and a lower KGE (0.57, (Figure 5(d)).This is because this course dataset comprises snow-depth measurements at significantly higher elevations than any other considered dataset (see elevation distribution in Figure 5(f)).
At those elevations, both interactions between the snowpack and topography complicate snow distribution compared to lowelevation areas, and density of input-data stations is much lower (Avanzi et al., 2021).Still, the fact that S3M does not show any obvious over-or underestimation for those elevations is encouraging for our scopes, also considering the comparatively coarse resolution of our implementation in Aosta valley (∼240 m).-water years).In this regard, Kling-Gupta Efficiency for snow depth and SWE is ∼ 0.70 and ∼0.84, respectively.Discrepancies between observed and simulated snow depth increase in spring, mainly because spatial heterogeneity in snow depth increases during the snowmelt season (Grünewald et al., 2010) and this challenges the comparison between a point measurement of snow depth and our 240-m snow model.Overall, the model correctly captures snowmelt rate and peak SWE for most water years, which are the primary variable of interest in operational snow hydrology.

Evaluation
Sporadic, yet abrupt oscillations in snow depth or SWE in the assimilated simulation are due to the assimilation dataset, which is the operational result of a multiregression model fitted across observed snow depth at ultrasonic-sensor stations and a number of physiographic features (Avanzi et al., 2021).These regressions often maintain -or even propagate -measurement noise, a frequent issue of ultrasonic snow-depth sensors (Ryan et al., 2008).On the other hand, open loop simulations do not display the same abrupt oscillations, which validates model's parameterizations.
S3M also reproduces the magnitude of bulk-snow density and its increase with time for all water years (Figure 6), in agreement with previous models implementing a similar parametrization of snow settling (De Michele et al., 2013;Avanzi et al., 2016) and despite its one-layer approach.Values of bulk and dry snow density are very close to each other during the accumulation season, while the latter diverges from the former during the snowmelt season.This is due to an increase in mass for the wet component of the snowpack during spring, as confirmed in terms of bulk liquid-water content (Figure 6).The seasonal range of variability for modeled bulk liquid-water content and its peak around 5 vol% during the snowmelt season agree with measurements by Techel and Pielmeier (2011); Heilig et al. (2015); Avanzi et al. (2017) and the international classification by Fierz et al. (2009).
The Torgnon study plot also measures incoming and reflected short-wave radiation, which allowed a comparison in terms of measured and simulated albedo (Figure 6).During the accumulation season, measured albedo is generally higher than simulated albedo; in particular, both measured and simulated albedo show maximum values around 0.95, but only the latter decreases well below 0.8-0.7 between snowfall events.Simulations by SNOWPACK at the same study plot (Terzago et al., 2019, not reported for brevity) qualitatively showed higher values than S3M, evidence that only relying on time as a predictor of albedo may yield frequent underestimations compared to a model that considers a broader spectrum of albedo predictors    Bulk-snow density was spatially fairly homogeneous, especially at the beginning of the snow season (Figure 7(f)).With time, some differences emerged, with snow density increasing faster in areas with both larger SWE and likely warmer temperatures (Figure 7(g) and (h)).Both the magnitude of snow-density values in Figure 7 and the fact that this variable was spatially more homogeneous than SWE tallies with previous works (López Moreno et al., 2013).

Evaluation: snow distribution
Maps of bulk-liquid water content were largely influenced by local climate, with a general rise in wetness with decreasing elevation that closely followed local topographic contours and aspect (Figure 7(i) to (k) and see Figure 8).Overall, bulk liquidwater content around the snow line was larger in April than in December or March, which again tallies with expected seasonal trends in wetness (Techel and Pielmeier, 2011).While liquid water in snow has been investigated for a long time (Colbeck, 1971), recent work on wet snow provides new opportunities for considering bulk liquid-water content from an operational standpoint, whether to predict the onset of snowpack runoff (Wever et al., 2014) or wet-snow avalanches (Mitterer et al., 2013;Wever et al., 2016).Early evidence that wet-snow conditions have increased in frequency and have extended well into the winter season due to a warming trend (Pielmeier et al., 2013) further justifies interest in wet snow predictions like those in Figure 8.
While application to significantly changed conditions will necessarily require some forms of parameter tuning or recalibration, S3M is among the few parsimonious snow models to provide spatially explicit, raster-like predictions of wet-snow patterns.
Maps of albedo showed the expected homogeneity during the accumulation season, owing to frequent snowfall events between November and March (Figure 7(l) and (m)).On the other hand, albedo in spring was much lower and spatially more diverse than in winter, with values larger than 0.8 at high elevation and values well below 0.6-0.7 in areas covered by older and wetter snow.

Evaluation: glacier evolution
To evaluate the glacier component of S3M in Aosta valley, we considered 94 ablation-stake measurements across the Rutor, Timorion, and Petit Grapillon glaciers (period 2009-2015, Figure 9).These are all high-elevation, debris free glaciers of various size (7.91 km 2 , 0.48 km 2 , and 0.18 km 2 , respectively -2012 data from the Aosta Valley Environmental Protection Agency), thus providing a representative sample to test the accuracy of S3M in capturing glacier ablation (note that these measurements were not included in the calibration of S3M, so they are a validation dataset).Simulations using the G3 glacier module (that is, including the ∆h parametrization) returned a correlation between simulated and observed change in thickness of 0.6 (Figure 9(d)).The correspondence between simulated and observed change in thickness across the range of variability in measurements was higher for Rutor than for Timorion and Petit Grapillon, which we interpret because of the large size of the former compared to model resolution (240 m).Also, the number of available samples for the Rutor is significantly larger than for the Timorion and the Petit Grapillon; for the latter, glacier surveys are challenged by extensive and deep crevasses, as well as frequent avalanches, that make point measurements non-representative of the actual melt pattern.Thus, the performance observed for the Rutor glacier is more representative of S3M predictive skills.
The evolution of glacier thickness for the Rutor glacier shows expected spatial patterns, with minor ablation at elevations above 3000 m ASL and progressively more intense melt close to the terminus below 2750 m ASL (Figure 10, see Figure 3 for location of this glacier in the study region).Annual changes show significant interannual variability, with somewhat 640 more intense melt as the evaluation period progresses; in any case, the spatial pattern is preserved as hypothesized by the ∆h parametrization.At elevations above 3200 m ASL, glacier thickness increased owing to conversion of seasonal snow to glacier ice at the end of the water year (see Section 2.6). Figure 11 shows similar spatial patterns for one of the most complex glaciers in the Alps, the Miage glacier, a 10.8-km 2 (2012), 10-km+ long valley glacier covering a ∼ 2000-m elevation range of the Mont Blanc massif (see Figure 3 for location of this glacier in the study region).Like many other valley glaciers across the southern Alps (Diolaiuti et al., 2003), vast portions of the Miage tongue are covered by debris, which has been shown to lead to below-debris melt being insensitive to variations in atmospheric temperature (Brock et al., 2010).Albeit hard to validate due to a lack of measurements, our implementation of the Miage glacier qualitatively captured this disconnection between intense melt across medium-elevation areas with little to none known debris and low melt rates in areas close to the glacier terminus that are well known to be covered by thick debris (Figure 11).Thanks to supporting a spatially explicit debris-driven melt factor (Equation 40), S3M yielded estimates of thickness change that were more spatially diverse and less correlated with elevation on the Miage than on the Rutor glacier (correlation coeffients of 0.95 and 0.85, respectively, compare Figure 11 with Figure 10).At elevations above ∼3000 m ASL, debris is residual if non-existent, and so change in thickness and elevation maintained the same high correlation observed on the Rutor glacier in Figure 10.Background layer is from the ESRI Satellite theme.

Applicability and future developments
As a hydrology-oriented cryospheric model, S3M delivers timely and computationally efficient predictions, while aiming at including the most salient processes of snow and glacier hydrology.Understanding this trade-off and its implications is important to determine model applicability and adequacy in a given context.
The structure and state variables of the model are all geared toward providing decision-relevant information for watersupply forecasting, which remains the prime area of application tested by the authors (Laiolo et al., 2014).In this context, typical questions that S3M contributes to answer are: how much snow-water resources are currently accumulated across the landscape?What headwater regions are currently releasing meltwater, and what are still accumulating snowpack?When is a certain percentage of the seasonal freshet expected to reach a given closure section?Are glaciers currently contributing runoff, and if so how much is their relative contribution?Decisions that are currently being informed by S3M thus range from flood-665 forecasting early warning to hydropower planning.S3M also provides some support to avalanche hazard forecasting, but this remains an unexplored area of application for which the microstructural detail in this model is largely insufficient.
Recently, we also started using S3M to produce future scenarios of water resources in mountain regions.Two particular features of S3M in this regard are the comparatively limited computational times of this model and the inclusion of both snow and glacier mass balances.Regarding the former, computational time for one water year worth of simulations in Aosta valley runs in ∼1.5 hours for a laptop with six i7 8th-generation cores and solid-state storage (240 m resolution, ∼3290 km 2 ).Note that S3M includes access to, and creation of, input and output NetCDF files (see the Supporting Information), and the frequency and size of output files in particular play an important role in determining computational time.Regarding snow and glacier mass balances, models providing even medium-level physical realism of both these features are still rare (Bongio et al., 2016;Li et al., 2015).Thanks to its integration with the Continuum hydrologic model (also open source, see Silvestro et al., 2013), S3M can deliver mass-conserving and spatially consistent predictions of the entire mountain water budget.
S3M is currently being actively maintained and further developed, with five main areas of planned future work.The first is the inclusion of wind effects, both as an additional component of the energy balance and as a driver of snow redistribution (wind drift).Explicitly including wind in the energy budget will help to resolve rain-on-snow events in the model, given that turbulent fluxes are a key contributor to melt during such flood-generating events (Marks et al., 1998;Würzer et al., 2016).As for wind drift, advances in this context are warranted not only in terms of relocation of blowing snow in the form of suspension, saltation, and creep (DeWalle and Rango, 2011), but also (and likely more importantly) in terms of the associated sublimation.
Progress in this regard has been hampered by a lack of detailed measurements of wind across the complex terrain of mountain headwaters, but recent datasets in this regard may favor future work on this topic (Guyomarc'h et al., 2019).
A second area of planned future work regards the inclusion of vegetation effects on the snowpack.The science of canopysnow interactions has identified four mechanisms through which vegetation can alter snowpack evolution compared to open areas: precipitation interception and throughfall, shortwave radiation shadowing, longwave-radiation enhancement, and wind shielding (Rutter et al., 2009).While the importance of each of these mechanisms for the fate of a seasonal snowpack dramatically changes with local climate (Lundquist et al., 2013), scientific consensus is that canopy may reduce peak SWE by more than 50% and lead to perturbations in the melt-out date on the order of weeks (Rutter et al., 2009), depending on canopy or snow-fractional cover.Helbig et al. (2020); Mazzotti et al. (2021) have recently proposed a parsimonious parametriza.tionfor canopy effects in large scale models, thus providing a solid starting point to include these processes in S3M.
The third direction of future development is liquid-water transport in snow, a rarely parametrized but important connection between surface melt and snowpack runoff.Water infiltration through snow manifests itself as both spatially homogeneous matrix flow and spatially heterogeneous preferential flow (Katsushima et al., 2013), with transition between these two regimes being driven by wet-snow metamorphism and snow properties like density and grain size (Avanzi et al., 2017;Hirashima et al., 2019).While capturing such micro-scale mechanisms is beyond the scope of a large-scale, distributed model, including some forms of preferential flow in S3M will likely enhance its performance in terms of timing and peak of early-season snowmelt events or rain on snow (Wever et al., 2014;Würzer et al., 2017).A way forward in this regard may be the simple parametrization originally proposed by Katsushima et al. (2009), which models preferential-flow discharge as a θ W -driven threshold process.
Another important aspect here is slope flow, that is, the tendency of snow to redistribute meltwater along layer boundaries for distances of hundreds of meters downhill (Eiriksson et al., 2013;Webb et al., 2018).Both measurements and consequently parametrizations of this process are still very rare, and more work is needed before this process can be included in parsimonious models like S3M.
Fourth, the conversion from snow to ice in S3M is very simplified, and completely skips the intermediate stage of firn (Cuffey and Paterson, 2010).While we usually turn off the glacier-mass balance when using S3M in flood forecasting, and while accumulation of multi-year snow as firn will likely characterize only very high elevations in a warming climate, this transition is still an important process to capture from the perspective of physical realism.Considering firn may also extend the applicability of S3M to polar regions, where for example firn-storage capacity is an important factor in determining the long-term fate of the Greenland ice sheet (Forster et al., 2014;Machguth et al., 2016).In this regard, Banfi and De Michele (2021) have recently proposed a local model of snow-firn transition for a binary-mixture snowpack like that considered in S3M.
Fifth, more work is planned to improve the data-assimilation component of S3M, particularly in terms of widening the array of state variables for which assimilation is supported (currently, SWE and snow depth) and in terms of improving assimilation data themselves.Regarding the former, the field of multivariate data assimilation is ripe (Piazzi et al., 2018), although scaling up these approaches across the landscape may come with significant computational requirements.Regarding new assimilation data, statistical learning is making promising steps towards mining new information from traditional, sometimes even sparse and fuzzy data across geosciences (Avanzi et al., 2019;Ghanjkhanlo et al., 2020;Revuelto et al., 2020;Dramsch, 2020;Grossi et al., 2021;Maurer et al., 2021;Shen et al., 2021;Mosaffa et al., 2022).Meanwhile, Cluzet et al. (2020) has showed some degree of success in assimilating satellite reflectance into snowpack simulations as a way to better capture snow microstructure and so the energy balance.Lessons learned from these recent advances will be incorporated in future releases of S3M.

Conclusions
We presented S3M v5.1, a spatially explicit hydrology-oriented cryospheric model that successfully reconstructs seasonal snow and glacier evolution through time.The model comprises parametrizations for precipitation-phase partitioning, snow and glacier mass balances, snow rheology and density evolution, snow aging and albedo, a hybrid temperature-index and radiationdriven melt parametrization, and provisions for data assimilation.Through a case study in Aosta valley, we showed that S3M provides robust performance across nearly twenty water years, with no systematic over-or underestimation of snow depth and a satisfactory simulation of the timing of peak accumulation and so the onset of the snowmelt season, including the transition from dry to wet snow (Root Mean Square Errors between S3M-simulated snow depth and in-situ measurements between ∼12 and ∼120 cm depending on elevation, with Kling-Gupta Efficiencies between 0.55 and 0.82).Overall, the model channels elements from the state of the art in cryospheric sciences into a parsimonious and computationally efficient model.Regarding snow, specific elements of relative novelty in this regard are an explicit representation of snow liquid-water content and a hybrid physics-based and temperature-index approach to snowmelt that decouples the radiation-and temperature-driven contributions.
The glacier component also includes the well-known ∆h parametrization and the possibility to feed the model with a distributed debris-driven melt factor, both comparatively new approaches in the field.S3M provides an open-source platform to simulate snow and glacier dynamics with the necessary physical realism for hydrologic purposes.Together with the hydrologic model with M D and M W the mass of the dry (snow grains) and wet (interstitial liquid water) constituents, respectively.The mass of air is assumed to be negligible compared to M D and M W (De Michele et al., 2013).

Figure 1 .
Figure 1.Prognostic variables and main mass fluxes of S3M (see Section 2 for details).αS is snow albedo (-), S f and R f are snowfall and rainfall rate (mm ∆t −1 ), respectively, SW ED (mm) and ρD (kg m −3 ) are dry Snow Water Equivalent and dry bulk snow density, respectively, SW EW is wet Snow Water Equivalent (mm), R, M , and O are snow refreezing, melt, and outflow, respectively (mm ∆t −1 ),

Figure 2 .
Figure 2. Values of the modulating parameter converting shortwave radiation into actual melt (m rad , panel a) and the degree-day melt parameter (mr, panel b) as a function of mean air-temperature over the previous 10 days ( T10d ) and various values of two tuning parameters (m rad and m r , see figure legends for values).

Figure 3 .
Figure 3. Aosta valley in north-western Italy.Panel (a): topography and glaciology of this region, with location of all snow-evaluation data used in this paper; panels (b) and (c): zoom on two intensive measurement regions, the Rutor glacier and the headwaters of Valpelline catchment, respectively."HS sensors" are continuous-time snow-depth ultrasonic sensors, "Torgnon" is an intensive study plot with a variety of snow and weather datasets (see Section 3.3 and Terzago et al., 2019), "Aval.probes" denotes locations of snow-depth measurements for avalanche-forecasting purposes, "Weekly samples" denotes locations of ∼weekly snow-depth measurements collected mainly for watersupply forecasting, while "Peak-HS courses" are snow-depth measurements collected along transects of several kilometers for hydropowerforecasting purposes (see Avanzi et al., 2021, for more details).Panel (a) includes location of the Miage and Rutor glaciers, for which detailed evaluation results are reported in Section 3.5.

Figure 4 .
Figure 4. Calibration objective metrics as a function of parameter values, with dModFactorRadS being m rad and a1dArctUp being mr (see the Supporting Information for details on model's notation in the source code).Parameter m rad is adimensional and mr is in mm °C−1 day −1 .For each parameter value Θ, the objective metric was computed using both all data from all years in one sample and separately for each water year (all-years and yearly values, respectively).Q1, Q2, and Q3 are the first, second, and third quartiles of objective metrics across all water years.WY is water year.

Figure 5 .
Figure 5. Performance of S3M in simulating point snow depth, as measured by continuous-time snow-depth sensors (a), the snow-depth sensor at the intensive study plot of Torgnon, and manual snow-depth measurements taken for avalanche, hydropower, and water-supply forecasting (c to e).Because simulations were carried out in assimilation mode, and assimilated maps indirectly involved snow-depth sensors and weekly snow samples, performances for these datasets are reported only for reference.Panel (f) shows the elevation distribution of each dataset and of Aosta valley for context.RMSE is Root Mean Square Error, KGE is the Kling-Gupta Efficiency (Kling et al., 2012), and r is Pearson's correlation coefficient.The spatial resolution of snow courses (panel d) is much finer than that of the model (say, ∼60 m vs. ∼220 m, respectively), thus a number of snow-course data were compared to the same modeling value.This explains the occurrence of sharp, horizontal lines in panel (d).
: the Torgnon study plot S3M reproduces the timing of peak accumulation and so the onset of the snowmelt season at the Torgnon study plot, both in terms of snow depth (h S ) and SWE (Figure 6, this Figure includes a mixture of calibration -2013-2019 -and evaluation -2020

Figure 6 .
Figure 6.Comparison between measurements and modeling outputs at the Torgnon study plot.First, second, third, and fourth columns are snow depth (HS), SWE, bulk-snow density, and albedo, respectively.The third and fourth columns also report dry-snow density and bulk volumetric liquid-water content, respectively.Each row is one water year (2013 to 2020).OL and Assim.stand for open-loop and assimilated simulations, respectively.

Figure 7 .
Figure 7. Simulated reanalysis of the 2020 snow season: spatially average air temperature and total precipitation (a), spatially average SWE and bulk liquid-water content θW (b), and distribution of SWE (c-e), bulk snow density ρS (f-h), bulk liquid water content (i-k), and albedo αS (l-n) for three example dates: 15 December 2019 (first column), 6 March 2020 (middle column), and 20 April 2020 (third column).Statistics in panel (a) and (b) are spatially averages across the entire Aosta Valley region.

Figure 8 .
Figure 8. Map of bulk volumetric liquid-water content for Aosta valley (LWC, same as θW ), 20 April 2020 at 12PM. Background layer is from the ESRI Satellite theme.

Figure 7
Figure 7 shows a simulated reanalysis of the 2019-20 snow season, which exemplifies the information and level of details provided by S3M to forecasters (note that this snow season is part of the validation pool).The 2020 snow season started by the end of October 2019 (Figure 7(a)), with largely uninterrupted precipitation events between November and December 2019 leading to nearly 75% of spatially average SWE across Aosta valley being accumulated before January 2020 (Figure 7(b)).January and February were relatively dry months, with only one storm in mid-February and then another one between February and early March 2020.The snowmelt season started in April, even though ∼40% of spatially averaged SWE at high elevations persisted by the end of May 2020.The season was characterized by an alternation between cold and warm spells, which led to frequent melt-freeze cycles in the spatially averaged snowpack (Figure 7(b)).

Figure 9 .
Figure 9. Spatial distribution of glacier-ablation measurements across the Rutor (a), the Timorion (b), and the Petit Grapillon (c) glaciers in Aosta valley.Panel (d) shows a comparison between measured and simulated ablation (positive values mean a decrease in local ice thickness).The number of points in panel (d) does not correspond to the number of stake locations in panels (a) to (c) because of repeated measurements taken cross multiple water years as the same stake location.Background layer is from the ESRI Satellite theme.

Figure 10 .
Figure 10.Rutor glacier: spatial patterns of annual change in glacier thickness (left), cumulative change in glacier thickness between September 2003 and 2019 as a function of pixel elevation (right, top), and spatial distribution of this cumulative change (right, bottom).Background layer is from the ESRI Satellite theme.

Figure 11 .
Figure 11.Miage glacier: spatial patterns of annual change in glacier thickness (left), cumulative change in glacier thickness between September 2003 and 2019 as a function of pixel elevation (right, top), and spatial distribution of this cumulative change (right, bottom).