Articles | Volume 15, issue 2
Data description paper
08 Feb 2023
Data description paper |  | 08 Feb 2023

IT-SNOW: a snow reanalysis for Italy blending modeling, in situ data, and satellite observations (2010–2021)

Francesco Avanzi, Simone Gabellani, Fabio Delogu, Francesco Silvestro, Flavio Pignone, Giulia Bruno, Luca Pulvirenti, Giuseppe Squicciarino, Elisabetta Fiori, Lauro Rossi, Silvia Puca, Alexander Toniazzo, Pietro Giordano, Marco Falzacappa, Sara Ratto, Hervè Stevenin, Antonio Cardillo, Matteo Fioletti, Orietta Cazzuli, Edoardo Cremonese, Umberto Morra di Cella, and Luca Ferraris

We present IT-SNOW, a serially complete and multi-year snow reanalysis for Italy ( 301 × 103km2) – a transitional continental-to-Mediterranean region where snow plays an important but still poorly constrained societal and ecological role. IT-SNOW provides  500 m daily maps of snow water equivalent (SWE), snow depth, bulk snow density, and liquid water content for the initial period 1 September 2010–31 August 2021, with future updates envisaged on a regular basis. As the output of an operational chain employed in real-world civil protection applications (S3M Italy), IT-SNOW ingests input data from thousands of automatic weather stations, snow-covered-area maps from Sentinel-2, MODIS (Moderate Resolution Imaging Spectroradiometer), and H SAF products, as well as maps of snow depth from the spatialization of over 350 on-the-ground snow depth sensors. Validation using Sentinel-1-based maps of snow depth and a variety of independent, in situ snow data from three focus regions (Aosta Valley, Lombardy, and Molise) show little to no mean bias compared to the former, and root mean square errors are of the typical order of 30–60 cm and 90–300 mm for in situ, measured snow depth and snow water equivalent, respectively. Estimates of peak SWE by IT-SNOW are also well correlated with annual streamflow at the closure section of 102 basins across Italy (0.87), with ratios between peak water volume in snow and annual streamflow that are in line with expectations for this mixed rain–snow region (22 % on average and 12 % median). Examples of use allowed us to estimate 13.70 ± 4.9 Gm3 of water volume stored in snow across the Italian landscape at peak accumulation, which on average occurs on 4 March ± 10 d. Nearly 52 % of the mean seasonal SWE is accumulated across the Po river basin, followed by the Adige river (23 %), and central Apennines (5 %). IT-SNOW is freely available at (Avanzi et al.2022b) and can contribute to better constraining the role of snow for seasonal to annual water resources – a crucial endeavor in a warming and drier climate.

1 Introduction

The seasonal snow cover is a key modulator of global climate (Flanner et al.2011) and a primary source of freshwater for more than one-sixth of the world's population (Barnett et al.2005; Immerzeel et al.2020). Snow water resources play a particularly important role in Mediterranean, summer-dry regions, where winter accumulation and the following summer freshet provide highly needed runoff to support societies and ecosystems as their demand peaks and precipitation declines (Zanotti et al.2004; Bales et al.2006; Viviroli et al.2007). Snow-dominated regions include, among others, the U.S. mountainous west, where snow provides 53 % of total runoff (Li et al.2017), central Asia (Immerzeel et al.2010), the Andes (Soruco et al.2015), and the European Alps (Viviroli et al.2007), which are regions from where the benefits of snow propagate downstream and across the globe (Sturm et al.2017). The critical role of snow for water resources, energy management, and ecosystem services is at the foundation of one of the most recurring, simplest, and yet most elusive questions in mountain hydrology (Bales et al.2006; Margulis et al.2015; Sturm et al.2017): how much snow is accumulated across the landscape at any given time?

Despite major advances since the seminal 1906 field campaigns by James E. Church on Mount Rose (NV, USA), this quest to quantify snow amount and distribution remains wide open (Dozier et al.2016). In situ measurements from ultrasonic snow depth sensors (Ryan et al.2008) or snow pillows (Cox et al.1978) are only representative of point conditions, with extrapolation at larger scales being hindered by the striking spatial heterogeneity of the snowpack (Grünewald et al.2010; Grünewald and Lehning2015; De Michele et al.2016) and possible perturbations of in situ instrumentation to snow natural conditions (Malek et al.2017). Extrapolation may be assisted by measuring snow amount along courses (Rice and Bales2010) or at strategically chosen locations that are representative of large-scale patterns (Zhang et al.2017b), but these solutions still imply intense labor and a comparatively high budget. Remote sensing, whether in the form of airborne lidar (Kirchner et al.2014; Painter et al.2016), remotely piloted aircraft (Bühler et al.2016; De Michele et al.2016; Harder et al.2016; Avanzi et al.2018), or optical and microwave satellites (Dietz et al.2012; Gascoin et al.2019b), has recently gained popularity in this context, particularly because it allows one to capture the full spatial distribution of the snowpack (Blöschl1999; Lievens et al.2019). However, remote sensing techniques are limited by either comparatively long revisit times, small areal coverage, uncertainties related to complex morphology, high maintenance costs, or cloud coverage. Finally, snowpack distribution models can simulate snow amount at virtually any resolution, but uncertainties in input data and in process representations make estimates based solely on the modeling of limited value in operational snow hydrology (Tang and Lettenmaier2010; Pagano et al.2014; Avanzi et al.2020).

Reanalyses obtained by assimilating in situ and remote sensing data into models are progressively becoming the most frequent, and arguably the most successful, solution to estimate snow water resources. Recent examples of such reanalyses for snow are the snow water equivalent (SWE) product by Margulis et al. (2016) across Sierra Nevada, California, the hyper-resolution ensemble-based reanalysis applied in Switzerland by Fiddes et al. (2019), the meteorological and snow reanalysis across the French mountains by Vernay et al. (2022), the High Mountain Asia UCLA Daily Snow Reanalysis by Liu et al. (2021), or the Austrian reanalysis product by Olefs et al. (2020). Estimates of snow coverage and amount are also available through Earth system reanalyses like the ERA suite by ECMWF (, last access: 19 July 2022), the NASA Global Land Data Assimilation System (GLDAS; see, last access: 19 July 2022), or the Japanese 55-year Reanalysis (, last access: 19 July 2022), among many others. Despite inheriting some of the original uncertainty in the data and models, reanalysis products optimally combine data and models in reconciled estimates and provide consistent coverage in space and time, thus paving the way for a new generation of snow science.

We present IT-SNOW, a  500 m snow reanalysis providing estimates of snow patterns across Italy ( 301 × 103km2) – a topographically and climatically complex region including some of the highest peaks in Europe (the Alps and the Apennines) and partially snow-dominated, socio-economically relevant regions like the Po river basin or the central Apennines. To our knowledge, this is the first open, sub-kilometric, serially complete, and multi-year snow reanalysis providing information on snow depth and mass for the Italian territory. Thus, IT-SNOW fills an important scale gap between in situ measurements and climate models or satellite-based datasets at kilometric resolution – such as the already mentioned ERA suite at 9 km, the H SAF suite (, last access: 19 August 2022), the Twentieth Century Reanalysis Project at  200 km or more (Compo et al.2011), the National Center for Atmospheric Research (NCAR) Climate Forecast System Reanalysis at  50 km or more (Saha et al.2014), the NASA Modern-Era Retrospective analysis for Research and Applications (MERRA) reanalysis product at  50 km or more (Gelaro et al.2017), or the non-mountainous Globsnow product at 25 km (Pulliainen et al.2020, see a complete review at, last access: 19 August 2022).

Figure 1Schematic of methods and data flows in S3M Italy, the operational chain used to generate the IT-SNOW dataset. All components of this chain are available in an open-source framework at (last access: 30 August 2022). SWE is snow water equivalent, LWC is liquid water content of snow, and SCA is snow-covered area.

IT-SNOW blends modeling, in situ data from snow depth sensors, and satellite observations from Sentinel-2, MODIS, and the H SAF initiatives and is the output of a real-time operational monitoring chain developed and maintained by CIMA Research Foundation for the Italian Civil Protection Department called S3M Italy (S3M stands for Snow Multidata Mapping and Modeling and is the underlying model used in this operational chain). The present dataset (IT-SNOW v1.0) includes daily reanalyzed outputs of SWE, snow depth, density, and bulk liquid water content from S3M Italy for water years 2011 through 2021 (a water year is defined as a period between 1 September and the following 31 August and is indicated with the calendar year in which it ends). Future updates are planned to expand this dataset (see Sect. 4). IT-SNOW is freely available at (Avanzi et al.2022b).

The paper is organized as follows. Section 2 describes S3M Italy (Sect. 2.1) and the preparation of the IT-SNOW reanalysis over the historical period 1 September 2010–31 August 2021 (Sect. 2.2). Section 3 evaluates the performance of IT-SNOW by using remote sensing data, in situ data, and an indirect water balance approach using streamflow. This section also includes a discussion of IT-SNOW sources of uncertainty (Sect. 3.3). Finally, Sect. 4 provides examples of use, while Sect. 5 details the data format and standards.

2 S3M Italy and IT-SNOW

2.1 The S3M Italy operational chain

S3M Italy provides real-time, spatially explicit estimates of snow cover patterns at  200 m resolution and with a latency of a few hours for the whole of the Italian territory ( 301 × 103km2). This operational chain includes algorithms to ingest in situ weather station data and satellite maps of snow cover, spatialization, and remapping tools to generate weather input and assimilation maps, parallel scripts to manage model simulations on multi-core servers, and a variety of post-processing and maintenance tools to generate final visualizations. S3M Italy is open source and available at (last access 30 August 2022), in particular through the Python package called fp-s3m. A schematic of the methods and data flows of this operational chain is reported in Fig. 1.

The core of S3M Italy (and hence of the reanalysis IT-SNOW) is the Snow Multidata Mapping and Modeling system or S3M (Avanzi et al.2022a). S3M is a spatially distributed cryospheric model solving the snow mass balance and parametrizing snowmelt using a hybrid temperature index and radiation-driven melt approach. Other processes included in S3M are snow settling, liquid water outflow, snow albedo evolution, and precipitation-phase partitioning. Complementary to these processes is an estimate of glacier melt based on the same hybrid approach used for snow but with modified parameters. Land cover effects, turbulent fluxes, and snow–forest interactions are currently not taken into account.

Figure 2Examples of input data used by S3M Italy to produce the IT-SNOW reanalysis. (a) Location of the radiation and relative humidity sensor stations for 4 February 2021 at 10:00 UTC. (b) Air temperature map for 14 February 2022 at 20:00 UTC, along with delineation of meteorological homogeneous regions in black (2019 update). (c) Precipitation map based on a modified conditional merging between precipitation gauges and radars for 23 November 2019 at 09:00 UTC (Bruno et al.2021). Background map from the Esri satellite theme. Note that some stations in panel (a) may host both types of measurements.

S3M is a raster-based model, where snow model equations are solved for each cell with no exchange of mass or energy across pixels. S3M is also open source and freely available at (last access: 30 August 2022), while more details on model physics and user requirements can be found in Avanzi et al. (2022a).

2.1.1 Input data preparation

Every hour at HH:40 UTC (where “HH” means every hour), the input data required by the model are downloaded and saved in predefined formats. These inputs include total precipitation, air temperature, relative humidity, and solar radiation, all of which are obtained from the database of the Italian Regional Administrations, Autonomous Provinces, and the Italian Civil Protection. Input data have an hourly time step. To fill potential gaps due to occasional malfunctioning and/or failures, every hour the automatic procedures check the existence of hourly inputs for the last 30 h. An unique estimate of the precision of the considered weather data is not available, as the type of sensor installed varies from one region to another. The installation and maintenance of the sensors generally follow guidelines from the World Meteorological Organization, to which the reader is referred (WMO2018).

Figure 3Daily national lapse rate climatology (a) and frequency distribution of daily national lapse rates according to S3M Italy for the period September 2010–August 2021. Q1, Q2, and Q3 are the first, second, and third quartiles of the daily national lapse rates for each month.


Total precipitation fields are the result of a modified conditional merging approach applied to precipitation gauges (spatial density of 1/100km2) and radar observations (Bruno et al.2021), so no further spatialization is performed in S3M Italy (Fig. 2). This modified conditional merging spatializes in situ precipitation data using an approach similar to kriging (called GRISO, from the Italian version of Random Generator of Spatial Interpolation from uncertain Observations), where the covariance structure is estimated for each precipitation gauge and each hour using radar data (see full details in  Sinclair and Pegram2005; Apicella et al.2021; Bruno et al.2021; Lagasio et al.2022). Final maps have a resolution of  1 km2, with a median root mean square error of less than 1 mm for a selection of 70 heavy precipitation events in Italy with accumulation greater than 100 mm or a maximum precipitation rate greater than 50 mm h−1 during the 2011–2014 period (see details in Bruno et al.2021). No phase partitioning is performed at this stage; separation between rainfall and snowfall is performed by the S3M model using the parametric approach by Froidurot et al. (2014), which relies on both air temperature and relative humidity.

Data of air temperature, solar radiation, and relative humidity are obtained as in situ point stations and further spatialized between HH:50 and HH+1:10 UTC ( 1 km for temperature and  500 m for radiation and relative humidity). For air temperature, the spatialization is performed by organizing station data into meteorological homogeneous regions as dictated by the Italian Civil Protection (see an example in Fig. 2; 2019 update) and fitting region-specific hourly linear regressions between air temperature and elevation (Figs. 2 and 3). These linear regressions are then applied using the meteorological homogeneous region's digital elevation model to derive temperature maps. Figure 3 reports monthly quartiles, and the frequency distribution of daily average national lapse rates, as derived through this procedure, which agree with the estimates by Rolland (2003) in the Alps. As for relative humidity and incoming shortwave radiation, we currently employ a computationally efficient method based on inverse distance weighting; no shadow effect or reflections from surrounding terrain are currently considered here, unless these are already captured by the comparatively dense network of stations.

Input data preparation ends between HH+1:10 and  HH+1:20 UTC, when input maps are cropped over the 20 computational domains, each corresponding to one Italian administrative region. Computational grids for these 20 domains were originally derived from a 20 m digital elevation model provided by the Italian Institute for Environmental Protection and Research (ISPRA), which was resampled at 200 m resolution using an averaging method. Besides elevation, S3M Italy employs static glacier maps from the Randolph Glacier Inventory v 6.0 (Pfeffer et al.2014).

Figure 4Data assimilated in S3M Italy to produce the IT-SNOW reanalysis. (a) Homogeneous snow regions based on snow climatology and expert knowledge. (b) Blended snow-covered area maps based on Sentinel-2, MODIS, and H SAF H10 snow products for 13 February 2022. (c) Snow-depth-interpolated map for 13 February 2022. For visualization reasons, the snow depth map in panel (c) was clipped using the concurrent snow-covered-area map. Snow depth estimates for some homogeneous regions were missing on that day due to insufficient in situ data; available regions are shown with 50 % transparency. Background map from the Esri satellite theme.

2.1.2 Assimilation data preparation

Data assimilation in S3M Italy is performed in the form of both satellite snow-covered area (SCA) and snow depth maps (Fig. 4). Maps of the snow-covered area are produced, once per day, by blending images from the ESA Sentinel-2, NASA MODIS, and EUMETSAT H SAF initiatives (product H10). First, 20 m maps of Sentinel-2 for the last 6 d across Italy are mosaicked (using the most recent one in case of overlapping images). A latency of 6 d is allowed both to manage the fairly infrequent revisit time of Sentinel-2 ( 5 d in Italy) and as a first provision to manage cloud obstruction. These maps at 20 m are then resampled at the  200 m grid of S3M Italy using the Python raster processing package Rasterio, with a mode resampling approach, to assign the dominant land cover among the 20 m pixels inscribed in each 200 m pixel. Once this first-guess SCA map is available, cloud-covered or unclassified pixels are further filled using resampled MODIS (, last access: 13 December 2022) and H SAF H10 data (, last access 13 December 2022), which have both nominal daily frequency and are used as distributed by the respective providers with no further processing. The result is a blended snow map providing information on snow cover, bare ground, and non-classified pixels. Besides mosaicking maps from multiple sources with different revisit times, no additional gap-filling for cloud coverage is performed. Assimilation for cloud-covered pixels was thus foregone. Like input data, this map is then remapped across the 20 regional domains to be assimilated.

Currently, Sentinel-2 SCA is produced by operationally applying the Sen2Cor algorithm by ESA (, last access 30 August 2022). Albeit lower in accuracy than snow-specific and high-resolution products like Theia (Gascoin et al.2019a), SCA maps derived with Sen2Cor were validated against snow depth sensors at the national scale during the average, representative, 2020 snow season and showed typical accuracy scores of the order of 0.7 to 0.8, as expected (see Fig. S1 in the Supplement and Main-Knorn et al.2017).

Snow depth maps are produced based on the interpolation of snow depth sensor in situ data. Every day during winter, measurements of the  350 snow depth sensors across Italy at 10:00 UTC are downloaded and quality-checked using an automatic filtering approach based on seasonality, climatological thresholds on minimum and maximum snow depth, and a filter based on a 6 h moving window coefficient of variation to detect grass growth after snowmelt. The specific device used to automatically monitor snow depth varies across the country, with the majority of them using an ultrasonic principle with an accuracy of ± 1 cm (Ryan et al.2008).

The  350 snow depth in situ data are then organized into 10 homogeneous regions defined along the boundaries of the International Standardized Mountain Subdivision of the Alps (SOIUSA; see Valt et al.2018), based on a tradeoff between maximizing data availability for each region and complying with expected climatology (e.g., we differentiated between inner-Alpine valleys and coastal, maritime, and mountain ranges; see Fig. 4 for a delimitation of these regions). For each of these homogeneous regions, a separate multilinear regression model is fitted across observed snow depth at sensor locations (predictand), elevation, slope, and aspect (predictors, with slope and aspect retained only if statistically significant). By applying the resulting (daily) multilinear regressions using a digital elevation model of each homogeneous region, daily snow depth maps are created and then cropped across the 20 computational domains of S3M. This is only done if at least 10 observations are available in a given homogeneous region; otherwise, spatialization and thus assimilation for that homogeneous region is foregone. An evaluation of this multilinear regression model in Aosta Valley showed biases of the order of 10 cm compared to avalanche probes, while a comparison with Sentinel-1 snow depth data at the national scale showed typical biases of up to 5 cm and typical root mean square errors below 10 cm (see Sect. 3 for details on these evaluation data).

Along with maps of snow depth, the procedure also generates a Kernel map quantifying spatial uncertainty in the multilinear regression model based on the distance across snow depth sensors (Avanzi et al.2021a, 2022a). This kernel is employed to assimilate snow depth map via a spatially distributed Newtonian relaxation approach (also known as nudging), under the assumption that uncertainty in snow depth maps will be lower in areas with a denser network of snow depth sensors. When a snow depth map is available, the assimilation procedure computes pixel-wise differences between a priori SWE and snow-depth-map-based SWE (after the conversion of snow depth maps into SWE maps using modeled density); this difference is then added to a priori SWE via Kernel weighting.

SCA maps are not assimilated directly but are used to clip snow-free pixels in snow depth maps before assimilation in the S3M model (thus assigning no snow to instances where snow depth maps estimate no snow although SCA maps observed snow). Both positive and negative differences are assimilated, meaning that the assimilation may result in either a decrease or an increase in simulated SWE. To further cope with grass interference in snow depth sensor data, the assimilation of snow depth and SCA maps is only performed between December and April, once per day, conventionally at 10:00 UTC.

2.1.3 Model runs and postprocessing

Upon completion of the input data remapping component of the modeling chain, parallel runs of S3M are performed (a first batch is launched at HH+1:24 UTC and a second one at HH+1:34 UTC). The runtime depends on the size of each modeling domain, with all simulations being completely roughly by HH+2:00 UTC; hence, there is a latency of less than 2 h for all domains. A Python wrapper manages each run, with S3M being a compiled Fortran executable. Every day at 3:00 UTC, summaries of previous day's simulations are compiled by mosaicking each domain on a national grid and saving outputs for visualization on the Italian Civil Protection WebGIS maintained by CIMA Foundation called myDEWETRA.

2.2 IT-SNOW preparation

For the scopes of IT-SNOW, we replicated an operational run of S3M Italy over the historical period from 1 September 2010 to 31 August 2021 (with a first period from 1 September 2009 through 31 August 2010 used as spinup). Historically observed weather data were thus downloaded and spatialized as outlined above, while also downloading, processing, and spatializing both satellite SCA maps and snow depth maps. Note that Sentinel-2 data were used only from summer 2021 and thus the assimilated SCAs before that period are the result of MODIS and H SAF maps. The native resolution of this historical run was  200 m, in line with the operational chain of S3M Italy.

Outputs of this historical run were saved every 6 h (at 5:00, 11:00, 17:00, and 23:00 UTC), and those at 11:00 UTC were assumed to be representative snapshots of daily conditions. These outputs at 11:00 UTC were thus remapped from the native  200 m grid to a national, geographic grid at  500 m (WGS84; EPSG 4326; pixel size of 0.005057). No projection was performed to avoid accuracy and distortion issues related to such cartographic systems. Owing to the geographic grid, the actual pixel size in meters changes with latitude, from a minimum of  460 m to a maximum of  508 m. We remapped at  500 m, using a nearest-neighbor approach for computational efficiency reasons and as an intermediate tradeoff maximizing on the predictive confidence between the native resolution of S3M Italy and the coarser resolution of precipitation data at 1 km2. Remapped outputs include instantaneous snow water equivalent (SWE), snow depth, bulk snow density, and bulk liquid water content, which overall form the IT-SNOW reanalysis dataset (see Sect. 5). Additional outputs are available upon request, and remapping over alternative grids is also possible.

3 IT-SNOW evaluation

3.1 Methods

Given that snow depth sensor data from the 20 regional networks were assimilated in IT-SNOW, we looked for alternative data that could act as truly independent evaluation sources. The first one was the satellite product by Lievens et al. (2019), C-SNOW, which provides daily, 1 km snapshots of snow depth across mountainous areas of the Northern Hemisphere based on an empirical change detection method applied to the Sentinel-1 measurements of the cross-polarization ratio. While C-SNOW is a remote sensing product, and as such may have locally larger uncertainties than in situ data, it has been successfully compared with in situ data from  4000 sites, with biases within ± 0.1 m for most of them (Lievens et al.2019). For the scope of evaluating IT-SNOW, the main advantage of C-SNOW compared to in situ data is that it is natively spatially distributed and thus allowed us to compare IT-SNOW across the landscape rather than at specific points. The evaluation period went from 1 September 2016 to 8 April 2020, which is the full span of 1 km C-SNOW data that are currently available. We remapped C-SNOW data onto the  500 m grid of IT-SNOW using a nearest-neighbor approach and then computed the pixel-wise bias and root mean square error with regard to IT-SNOW. Note that C-SNOW data are available only for dry snow conditions and that the optimization procedure included some (but not all) of the Italian in situ snow depth sensor data.

The second source of validation data considered here were in situ data taken in Aosta Valley (northwestern Italy; see Fig. 6), Lombardy (northern Italy; see Fig. 7), and Molise (central Italy; see Fig. 8). These three areas present significantly different climates and thus snow types (Sturm and Liston2021), with Aosta Valley and Lombardy being characterized by a seasonally consistent and deep Alpine snowpack and Molise more exposed to lake effect snowfalls from the Adriatic Sea and thus to a more ephemeral and maritime snow cover (Da Ronco et al.2020). These three datasets are topographically diverse and cover a comparatively long time span (see below), meaning that they are representative of larger-scale performances of IT-SNOW.

Measurements in Aosta Valley included yearly snow course manual samples taken at peak accumulation every 50–100 m along elevation transects of several kilometers upstream of five hydropower reservoirs (water years 2011–2021; elevation above 2000 ma.s.l., above sea level), daily to weekly manual measurements at recurring locations for avalanche forecasting (water years 2011–2021; elevations above 1000 ma.s.l.), and hourly automatic measurements from an ultrasonic snow depth sensor and a SWE sensor in Torgnon (2012–2020; elevation 2160 ma.s.l.; see a data inventory in Avanzi et al.2021a, 2022a). Data in Lombardy include weekly estimates of snow depth and SWE obtained by running the physics-based multi-layer SNOWPACK model (Bartelt and Lehning2002) in correspondence with automatic weather and snow stations at medium elevation (1800–2600 ma.s.l.), where the model was forced using local input data and assimilating local snow depth (henceforth AWS snow depth and SWE – years 2016 through 2021, with AWS standing for automatic weather station) and measurements of snow depth and SWE collected between May and June on glacier terrains at a very high elevation for mass balance purposes (elevations above 3000 ma.s.l.; years 2016 through 2021; henceforth glacier snow depth and SWE). Data in Molise included daily to weekly manual measurements at four recurring locations for avalanche and water supply forecasting (2011–2021; elevation 1200–1500 ma.s.l.). Performance metrics between observed and simulated SWE, snow depth, and bulk snow density included bias, root mean square error, mean absolute error (MAE), positive and negative mean error (PME and NME, respectively), the Kling–Gupta efficiency (Kling et al.2012), and Pearson's correlation coefficient.

The third source of validation data were streamflow measurements for a selection of 102 basins in Italy for which long-term, serially complete, and quality-checked time series of streamflow were available for the period from 1 September 2010 through 31 August 2019 (Bruno et al.2022). We used these data to compare the annual peak of water stored in snow and annual cumulative streamflow at the closure section of these basins (both in Gm3) as a proxy for the proportion of annual flow that was accumulated as snow. To this end, water stored in snow (simply SWE in Gm3 in the following) was obtained by multiplying pixel-wise SWE (in mw.e.) from IT-SNOW by the area of each cell and then summing all pixel-wise values. Given that it is general knowledge that Italian precipitation climatology is a mix between snow and rain, we not only expect these ratios to be between 0 and 1 but also to predominately be smaller than 0.5. Owing to precipitation increasing with elevation (Avanzi et al.2021a), we also expect these ratios to increase with the average elevation of each of these catchments. While indirect in that IT-SNOW is not evaluated against snow data, this third evaluation stems from a long-standing tradition of inverting the hydrological cycle (Valery et al.2009) to provide insights into the consistency of IT-SNOW estimates with the local water budget.

Figure 5Evaluation of IT-SNOW using the Sentinel-1 snow depth product C-SNOW for the period September 2016 through April 2020. Pixel-wise bias for the Italian Alps (a) and the central Apennines (b), frequency distributions of pixel-wise bias (c) and root mean square error (d) are shown. The background map is from the Esri satellite theme. Panels (a) and (b) refer to the two areas of Italy with seasonally deep snow cover.

3.2 Results

3.2.1 Evaluation of snow depth against C-SNOW

Mean pixel-wise bias between IT-SNOW and C-SNOW was close to zero (0.01 m), with a median value of zero, and the first and third quartiles being 0.03 and +0.02 m, respectively (Fig. 5c). Thus, the distribution of spatial biases was well centered around zero (Fig. 5c). Mean pixel-wise RMSE was 0.22 m (Fig. 5d), i.e., close to the mean absolute error found by Lievens et al. (2019) in the original evaluation of C-SNOW with snow depth sensors (0.18 to 0.31 m). The first, second, and third quartiles of pixel-wise RMSE were 0.14, 0.19, and 0.27 m, with only a fraction of values above 0.5 m (Fig. 5d; both bias and RMSE were calculated between time series at each pixel). We conclude that the two products provide consistent estimates of snow depth across the Italian landscape.

The spatial distribution of bias across the Italian Alps and the central Apennines showed no obvious pattern, with only a tendency of IT-SNOW to underestimate C-SNOW snow depth at high elevations (see Fig. 5a and b). This may be related to well-known biases of precipitation gauges and radars, the main sources of input precipitation in S3M Italy, at high elevations and/or in inner-Alpine areas (Zhang et al.2017a; Cui et al.2020; Avanzi et al.2021a, see Sect. 3.3 for a discussion). Besides this bias with elevation, biases and RMSEs between IT-SNOW and C-SNOW were consistent across the 10 homogeneous regions used to generate snow depth maps to be assimilated in IT-SNOW (Figs. S2 and S3 in the Supplement).

Figure 6Evaluation of IT-SNOW in Aosta Valley. (a) Topography of the focus region, along with sampling location of evaluation data. Panels (b) and (c) show observed vs. simulated snow depth at snow course and avalanche probe locations. Panels (d) and (e) show simulated vs. observed snow depth at Torgnon. Panels (f) and (g) show simulated vs. observed SWE at Torgnon. Panels (h) and (i) show simulated vs. observed bulk snow density at Torgnon. Note that simulated snow depth, SWE, and density in panels (d), (f), and (h) were smoothed using a 5 d moving window for readability. r is Pearson's correlation coefficient, HS is snow depth, and SWE is snow water equivalent.

Table 1Overview of IT-SNOW performance versus in situ data in Aosta Valley, Lombardy, and Molise. RMSE is the root mean square error, MAE is the mean absolute error, PME is the positive mean error, NME is the negative mean error, KGE is the Kling–Gupta efficiency (Kling et al.2012), and r is Pearson's correlation coefficient.

Download Print Version | Download XLSX

Figure 7Evaluation of IT-SNOW in Lombardy. (a) Topography of the focus region, along with location of evaluation data and an example of glacier data sampling geometry (b, c) IT-SNOW vs. SNOWPACK snow depth and SWE at medium-elevation snow stations, respectively. (d, e) Observed vs. IT-SNOW snow depth and SWE at very-high-elevation glacier sites in the context of end-of-season mass balance surveys. r is Pearson's correlation coefficient, AWS is automatic weather station, and SWE is snow water equivalent.

3.2.2 Evaluation against in situ data

IT-SNOW estimates of snow depth, SWE, and bulk snow density were generally well correlated with measurements in Torgnon (Aosta Valley; Fig. 6e, g, and i), with local RMSE for snow depth, SWE, and density being 30 cm, 95 mm, and 93 kg m−3, respectively (Table 1), and Pearson's correlation between 0.43 and 0.81 (Table 1). Biases were minor in Torgnon, with Kling–Gupta efficiencies that were significantly above the no-skill threshold of 0.41 (see Knoben et al.2019, and Table 1). From a seasonal perspective, IT-SNOW and measurements in Torgnon also agreed in terms of accumulation and snowmelt temporal patterns, in addition to the date of peak accumulation (Fig. 6d, f, and h). The lowest correlation was found for bulk snow density (Table 1), which is not surprising, given that density was indirectly derived from measurements of snow depth and SWE (DeWalle and Rango2011), a procedure that may increase noise (Terzago et al.2020). Performance scores decreased when considering avalanche probes (Fig. 6c and Table 1) and high-elevation snow courses (Fig. 6b and Table 1), which was expected, given the significant scale of mismatch between a  500 m reanalysis and topographically diverse, in situ manual measurements and the already mentioned possible underestimation of precipitation fields at high elevations (see Fig. 5 and  Avanzi et al.2021a). Overall, evaluation results in Aosta Valley showed that IT-SNOW successfully reconstructs both the seasonal dynamics and peak timing of snow depth and mass, and hence density, of this heavily instrumented region (see Figs. 2 and 4).

In Lombardy, IT-SNOW estimates of snow depth and SWE generally agreed well with those provided by SNOWPACK at medium elevations, despite somewhat larger RMSEs and biases, larger mean absolute errors, and a lower correlation than in Aosta Valley (Table 1 and Fig. 7). Note that these medium-elevation data in Lombardy were not directly observed but were the result of modeling, which may have increased their own uncertainty. On the other hand, IT-SNOW in Lombardy showed an expected systematic underestimation of very-high-elevation mass balance measurements on glaciers (biases of 717 mm and 110 cm; Table 1), despite a promising correlation and a KGE between observations and simulations (Table 1). While overcoming these underestimations is a major challenge for large-scale reanalyses like IT-SNOW, again because of the scale and undercatch issues discussed with regard to C-SNOW (Fig. 5) and Aosta Valley data (Fig. 6), and while very-high-elevation regions play only a minor role in a catchment water balance, we discuss pathways to improve IT-SNOW in this regard in Sect. 3.3.

Figure 8Evaluation of IT-SNOW in Molise. (a) Topography of the focus region, along with sampling location of evaluation data. (b–d) Simulated vs. observed snow depth, SWE, and bulk snow density at Molise snow stations, respectively. (e–g) Example of simulated vs. observed time series of snow depth, SWE, and bulk snow density at Pescopennataro, respectively. Note that simulated snow depth, SWE, and density in panels (e–g) were smoothed using a 5 d moving window for readability. r is Pearson's correlation coefficient, HS is snow depth, and SWE is snow water equivalent. Density was measured with a resolution of 20 kg m−3; hence, the discrete values of the point cloud along the x axis are shown in panel (d).

Correlations between measured and simulated snow depth, SWE, and density were lower in Molise than in Aosta Valley and Lombardy (Fig. 8b, c, and d; RMSE of 61 cm, 200 mm, and 109 kg m−3, respectively, and biases of 34 cm, 82 mm, and 24 kg m−3, respectively; Table 1). We interpret this outcome as being due to the sparser network of assimilated snow depth sensors used in southern Italy compared to northern Italy (see Fig. 4). Yet, IT-SNOW successfully reconstructed not only the seasonal dynamics of accumulation of melt in Molise but also the much more significant interannual variability in peak SWE than in Aosta Valley (Fig. 8f).

These performances of IT-SNOW against in situ data are consistent with snow reanalyses over other areas of the Alps. In this regard, Vernay et al. (2022) report median RMSEs for snow depth of the order of 10–40 cm in France (maximum RMSEs up to 90 cm), with an increasing trend with elevation; both results echo findings in this paper for areas at medium elevation, where the bulk of forcing and evaluation data are available (Avanzi et al.2021a). In Austria, Olefs et al. (2020) found RMSEs for snow depth and SWE of the order of 10 cm and 100 mm, with correlations of 0.86 and 0.91, respectively; this accuracy is higher than that of IT-SNOW, likely because of the much more homogeneous coverage of snow data in Austria compared to Italy (see Fig. 1 in Olefs et al.2020). Finally, Fiddes et al. (2019) report RMSEs of the order of 38–53 cm for snow depth and 184–258 mm for SWE, which again tallies with the accuracy of IT-SNOW.

Figure 9Indirect validation of IT-SNOW with streamflow data. (a) Mean annual ratios (in %) between peak SWE and annual cumulative streamflow for a selection of 102 water basins with long-term, serially complete, and quality-checked time series of streamflow. (b) Frequency distribution of these annual ratios across all sections and years. (c) Correlation between annual streamflow and annual peak SWE across all sections and years (red line is the 1:1 reference line; the black dashed line is a linear regression between peak SWE and annual streamflow). (d) Correlation between mean peak SWE and river basin mean elevation (the black dashed line is a linear regression between these two data). Q is annual streamflow. Note that considered water sections of the central Po river valley do not account for snow accumulated in the Swiss canton Ticino, which is not included in IT-SNOW (Ticino represents about 5 % of the Po river basin at its most downstream closure section). The background map is from the Esri satellite theme. SWE is expressed here in the form of total volume of water stored in snow (Gm3; see Sect. 3.1).

Simulated time series of snow depth and SWE in Torgnon present frequent, abrupt oscillations that are not related to any physical process such as melt or settling (Fig. 6d and f). These oscillations, which were already noted in Avanzi et al. (2022a), are due to the assimilated snow depth maps, which often include – or even propagate – instrumental noise (Ryan et al.2008). Indeed, similar oscillations are not visible in bulk snow density time series, which are unaffected by assimilation (Fig. 6h), or in Molise (Fig. 8e and f), where assimilation is much rarer due to a sparse network of snow depth sensors. While these oscillations do not affect the seasonal reconstruction of snow depth or the temporal patterns of peak SWE, it is important to bear them in mind when performing temporally fine analyses with IT-SNOW. Such oscillations could be better handled by assimilation techniques that explicitly account for point measurement uncertainty, such as a Kalman or particle filter (Piazzi et al.2018), but doing so would entail significantly higher computational demands that are currently non-feasible given the tight, real-time schedule required by S3M Italy.

3.2.3 Comparison of estimated SWE with observed streamflow

On average, peak SWE is about 22 % ± 24 % of annual streamflow across the considered 102 gauge stations and 9 water years (median of 12.8 % for the period 2011–2019; see Fig. 9b. Note that SWE is expressed in this section in the form of total water volume stored in snow, Gm3; see Sect. 3.1). As expected, the distribution of this ratio is highly skewed toward values below 50 % and shows a clear latitudinal trend, with higher values for rivers draining from the southern Alps in the Po river valley and a handful of basins draining from central Apennines (Fig. 9a). In these snow-dominated regions, the average ratio across the considered 9 water years locally exceeds 50 % to 60 %, especially in the high-elevation regions of northwestern Italy and the Adige valley. Comparatively small watersheds on the central–western side of the Italian peninsula show significantly smaller values, reflecting lower elevations and a more maritime and warmer climate.

As an additional, indirect validation of IT-SNOW, peak SWE is significantly correlated with annual streamflow (Pearson's correlation coefficient of 0.87; Fig. 9c). In other words, peak SWE is a robust predictor of annual streamflow in Italy, which agrees with past findings in more snow-dominated regions of the world (Pagano et al.2004; Rosenberg et al.2011; Harrison and Bales2016). Again, as expected, peak SWE is smaller than annual streamflow (Fig. 9c) and is significantly correlated with the mean elevation of each watershed. Overall, these evaluations show consistency between IT-SNOW SWE estimates and the Italian water budget.

3.3 Sources of uncertainty

Like all reanalyses combining sparse data and a model over large domains, IT-SNOW is the result of a number of tradeoff choices and epistemic uncertainties that should be taken into consideration when handling this dataset. The first source of uncertainty is represented by precipitation estimates, similar to any snow–hydrologic model simulation in mountain regions. The modified conditional merging approach used in IT-SNOW has already been extensively validated and shows robust performances for heavy precipitation events at the national scale (relative error in high flows < 25 % for 72 % out of 241 Italian river sections when used to force a hydrologic model; see  Bruno et al.2021). However, it does not include explicit provisions for reconstructing precipitation orographic gradients (besides those captured by the location of stations and by radar images). Previous work in Aosta Valley (Avanzi et al.2021a) and elsewhere (Lundquist et al.2015; Zhang et al.2017a; Avanzi et al.2020), in addition to the results reported in this paper (Figs. 5 and 7), show that including these orographic gradients is important to close the water budget of small, high-elevation, Alpine catchments. Targeted validations for such high-elevation Alpine catchments will thus be the subject of future work.

A second source of uncertainty related to precipitation and more generally to in situ weather data is data sparsity and how this can affect spatialization in partially ungauged regions. While the density of weather data in Italy is comparatively high (as an order of magnitude, there is  1 station every 100 km2 or more), and stations are routinely monitored and maintained by regional administrative authorities, we expect estimates for regions at the boundary of the Italian territory to show an inherently larger uncertainty than the rest of the country (because of a lack of input data outside the Italian territory). Data sparsity also limits the amount of snow depth data that we can currently employ in assimilation and the extent of evaluation regions in this paper (see Sect. 3). While estimating a real-time layer of uncertainty for this reanalysis product is currently not feasible, this will be target of future work.

Regarding the assimilation data, uncertainty in satellite SCA is in line with standards by the European Space Agency (ESA), as already discussed in Sect. 2. On the other hand, we noted that noise in snow depth sensor data, along with the likely simplistic multilinear regression approach used to spatialize snow depth across the landscape, often introduces artifacts in snow depth that translate into abrupt oscillations in snow depth values at the daily timescale (Fig. 6). Another source of uncertainty in this regard is that we organized the Italian territory in 10 homogeneous regions, but no smoothing at the regional boundaries is currently in place. As a result, IT-SNOW may sometimes overestimate snow depth variability across the boundaries of these homogeneous regions and/or exaggerate the role of single topographic predictors of snow depth, such as aspect. An alternative to our approach would be to directly assimilate remote sensing products; ongoing efforts by the European Space Agency and others are moving towards this direction, and we aim to include new findings in this regard in our assimilation framework.

In terms of epistemic uncertainty, the S3M model relies on an enhanced temperature index approach that was calibrated in Aosta Valley and then extensively evaluated elsewhere in Italy (both in this paper and in other projects, such as Alfieri et al.2022). Relying on the same parameters across the whole country could introduce some additional uncertainty at the local scale. However, results in this paper show a credible reconstruction of melt dynamics, even in areas with only occasional assimilation (Fig. 8). In this regard, Bouamri et al. (2018) reported encouraging results when transferring model parameters of the same enhanced temperature index approach to uncalibrated sites. A more influencing factor in this sense could stem from S3M not solving the full energy balance, but Magnusson et al. (2015) show that this is not a key driver of model performance for snow bulk variables at the daily scale.

Figure 10Example of use of IT-SNOW. (a) Spatially aggregated total daily SWE across Italy. (b–l) Mean annual SWE (November to June for water years 2011 to 2021). The background map is from Esri terrain. SWE in panel (a) is expressed in the form of total volume of water stored in snow (Gm3; see Sect. 3.1).

4 Examples of use

By providing spatially explicit, high-resolution, and serially complete estimates of snow patterns, reanalyses like IT-SNOW have the potential to fill important knowledge gaps in hydrology, for example, by elucidating the role of snow in supporting worldwide water security (Viviroli et al.2007) or snow sensitivity to climate extremes like droughts (Hatchett and McEvoy2018). IT-SNOW also responds to pressing societal questions from operational water resources managers and decision-makers, who need diversified information on snow distribution, amount, and interannual variability to make accurate decisions regarding water allocation and use (Harrison and Bales2016). In this section, we show examples of how IT-SNOW can be used to answer the following two societally relevant questions (Dozier2011; Margulis et al.2015; Painter et al.2016): (1) how much snow is accumulated across the landscape? (2) Where is it distributed?

4.1 How much snow is accumulated across the landscape?

Figure 10a shows total daily SWE across Italy for the whole period of record, again in the form of total water volume in snow (Gm3; see Sect. 3.1). Peak SWE in Italy is on average 13.70 ± 4.9 Gm3, with a minimum in 2017 ( 5 Gm3) and a maximum in 2014 ( 23 Gm3). SWE peaks are, on average, on 4 March ± 10 d, i.e., about 1 month earlier than the 1 April reference date; this finding is in agreement with a recent reconsideration of this conventional date (Montoya et al.2014). The earliest peak SWE date was 3 February (2019), while the latest was 26 March in 2013. Seasonal dynamics show signs of intraseasonal melt (e.g., see late 2011 to early 2012), which is expected in a Mediterranean region where cold-alpine and maritime snow types coexist (Sturm et al.1995; Sturm and Liston2021). Nonetheless, snow seasons are temporally continuous, with no episode of intraseasonal melt out. Little to no carryover occurs between seasons.

A large fraction of Italian snow accumulates across the southern Alps, while snow accumulation on the Apennines is spatially more limited and – importantly – more variable from one season to the next (Fig. 10b to l). This increased variability in the Apennines compared to the Alps agrees with sparse, but consistent, previous work showing that snow in the Apennines – and particularly on their eastern side – is the result of intense, highly seasonal, lake effect storms (Da Ronco et al.2020). One such event is evident in the 2012 reanalysis of IT-SNOW (Fig. 10c), when central Italy was hit by extensive snowfalls as part of a continental cold wave (Demirtaş2017). Data for February 2012 report 150 cm or more of fresh snow in places, with more precise estimates of the return period for this event being challenged by the sparsity of the data network (Bisci et al.2012). Another similarly exceptional event was the 2014 season in the southern Alps (Fig. 10e), with 8 m of peak snow depth at 2000 m (200 % of the long-term mean; see Chiambretti et al.2014). IT-SNOW also correctly captured rare but significant low-elevation snowfall events (see water year 2011; Fig. 10b).

The spatially and temporally consistent framework of IT-SNOW can aid not only climatological and hydrologic studies, for example, by providing highly necessary validation datasets of snow cover and SWE, but also operational water resources managers who routinely use SWE as an estimate of freshet volume (Harrison and Bales2016). Similar efforts have already been made in other regions like California, USA, where Margulis et al. (2016) estimate about 20 Gm3 of mean peak water volume in snow across the Sierra Nevada, or Japan, where Niwano et al. (2022) estimate 42.2 Gm3 on average during the 2017–2022 winters. These estimates tally with those of IT-SNOW across the Italian mountain ranges (Fig. 10), an area with a similar geographic span to the Sierra Nevada, California, or Japan but a less snow-dominated climate. With its high spatial resolution and fine temporal granularity, IT-SNOW can contribute to the quest for constraining the extent and volume of the world's cryosphere.

Figure 11Example of use of IT-SNOW. (a, b) Percentage of mean winter SWE that is accumulated across the major Italian river basins (right) and location of these major basins across Italy (left). (c–i) Median and quartiles of daily basin-wide SWE across a selection of the most snow-dominated basins in Italy. Q2 is the median, and Q1 and Q3 are the first and third quartiles, respectively. The background map is from the Esri satellite theme. SWE is expressed here in the form of total volume of water stored in snow (Gm3; see Sect. 3.1).

4.2 Where is it distributed?

Spatially aggregating mean winter SWE across major Italian watersheds shows that  52 % of Italian snow water resources accumulate across the Po river basin, the largest Italian water basin (Fig. 11b). The second-largest snow reservoir in Italy is, as expected, the Adige river basin (23 %), followed by the Piave, Tagliamento, and Brenta basins (6 %, 3 %, and 3 %, respectively). Collectively, these five Alpine water basins host nearly 87 % of Italian snow. A second hotspot for snow accumulation is the central Apennines, with the Tiber river basin accumulating about 2 % of the national mean winter SWE. This and other three basins in central Italy (Aterno–Pescara, Garigliano, and Sangro) accumulate about 5 % of national SWE, with the remaining 8 %–9 % scattered across the remaining basins.

Median and quartiles of daily SWE for a selection of the most snow-dominated Italian basins show the typical seasonal dynamics of snow accumulation and melt, with the peak SWE date between water year days  180 and  200 (that is, between the end of February and mid-March; see Fig. 11c to i. SWE is expressed as total volume of water stored in snow, Gm3; see Sect. 3.1). As already noted in Fig. 10, the Apennines basins show much more interannual variability than Alpine basins, as evident from the larger spread between the first and third quartiles in Fig. 11g to i compared to Fig. 11c to f.

Italian basins also show significantly different melt rates, e.g., the mean annual SWE for the Po river basin is 4 Gm3, with an almost continuous snow cover and an almost symmetrical accumulation and melt season (Fig. 11c). This symmetry significantly differs from more radiation-driven and maritime snow types, such as that of the Aterno–Pescara or Tagliamento rivers (Fig. 11f and h), where the melt season is shorter than the accumulation season.

With its initial time span of 11 years, IT-SNOW-derived statistics like those from Fig. 11c to i can provide estimates of medium-term SWE variability and thus put accumulation at any given time (like those in Fig. 10) into a broader context. Such a broader context is becoming all the more important in a warming and drier climate, particularly given the growing emergence of snow droughts as a specific typology of droughts (Harpold et al.2017; Hatchett et al.2017; Huning and AghaKouchak2020). In this regard, medium-term climatological bands like those in Fig. 11 have already been used to contextualize the severe 2022 precipitation deficit in northern Italy by showing that this deficit translated into  40 % SWE compared to the 2009–2021 median (Toreti et al.2022). While the initial time span of IT-SNOW snow is too short for rigorous deficit estimates, the operational chain delivering this reanalysis will provide yearly updates for future seasons, while an extension to the near past (say, 2002–2009) is also in consideration. Besides deficit analysis, we expect the notion of snow accumulation and melt rates to aid in other contexts, such as ecology (Slatyer et al.2022) or climate change assessments (Musselmann et al.2017).

5 Data format

IT-SNOW is available in an open-access framework at (CC BY-NC 4.0; see Avanzi et al.2022b). Data are organized in monthly netCDF files, with zlib compression, each hosting a 3D matrix of daily maps for one variable of interest (SWE, snow depth, bulk snow density, and bulk liquid water content). Filename strings follow a consistent convention, including variable name and month (e.g., for SWE data of September 2010). Variable labels and units are as follows: SWE for snow water equivalent (mmw.e.), HS for snow depth (cm,  Fierz et al.2009), RhoS for bulk snow density (kg m−3), and Theta W for bulk liquid water content (vol%).

In addition to data, each netCDF file includes information regarding the reference system (variable crs, including well-known text strings for EPSG 4326), latitude and longitude matrices, and a time array. Compatibility with the NASA Panoply system (, last access: 23 August 2022) and with QGIS were both verified, including reprojection to a UTM (Universal Transverse Mercator) metric system. Each netCDF file also includes metadata identifying contact points and curators.

6 Code and data availability

IT-SNOW is available in an open-access framework at (CC BY-NC 4.0; see Avanzi et al.2022b).

Sources of the data used are reported in the paper and include the database of the Italian Regional Administrations and Autonomous Provinces, as accessible by CIMA Research Foundation through the Italian Civil Protection, the Aosta Valley Environmental Protection Agency (snow courses), the Aosta Valley Ufficio Neve e Valanghe (avalanche probing data), Meteomont (as available for the Molise Region; Molise snow data), the Lombardy Environmental Protection Agency (Lombardy data), and the C-SNOW initiative (Lievens et al.2019).

The S3M model is open-source software. Official releases are available at (Avanzi et al.2021b). The most important components of the S3M Italy operational chain are also open source (see, Avanzi et al.2021c).

7 Conclusions

We presented IT-SNOW, a spatially explicit and multi-year reanalysis of snow cover patterns across Italy at  500 m resolution. IT-SNOW is the reanalyzed output of S3M Italy, a cryospheric modeling chain operationally delivering spatial snapshots of snow water resources for civil protection applications. Through S3M Italy, IT-SNOW ingests input data from thousands of automatic weather stations across the Italian territory, while assimilating daily snow-covered-area maps from ESA Sentinel-2, NASA MODIS, and EUMETSAT H SAF products and multilinear regressions of on-the-ground snow depth data. Validation results show little to no mean bias compared to C-SNOW, a state-of-the-art retrieval of snow depth from Sentinel 1, root mean square errors of the order of 30–60 cm for in situ measured snow depth and 90 to 300 mm for in situ measured snow water equivalent, a strong (0.87) correlation between peak SWE and annual streamflow, and ratios between peak SWE and annual streamflow that are in line with expectations for this mixed rain–snow region (22 % on average). Examples of use showed how IT-SNOW can both fill fundamental knowledge gaps in snow hydrology and support real-world applications by answering recurring questions like “how much snow is accumulated across the Italian landscape?” (on average 13.70 ± 4.9 Gm3 peak SWE), or “where is it?” ( 52 % and 21 % across the Po and Adige river basins, respectively, with the remainder between less snow-dominated watersheds in northeastern Italy and the central Apennines).

IT-SNOW will be updated annually, and community engagement will be favored to maintain a high-resolution reanalysis of snow for Italy. Engagement will follow two lines of action. The first is institutional, as we are engaging with the Italian regional administrations to benchmark IT-SNOW against their long-standing systems providing local to regional estimates of SWE with methodologies that vary among the offices. The second is bottom-up, as we are designing tools to favor engagement through a GitHub discussion page hosted at (last access: 30 August 2022). The author team thus remain open to critical feedback from the user community on IT-SNOW accuracy, issues, and improvement opportunities.


The supplement related to this article is available online at:

Author contributions

FA, SG, and FD developed S3M Italy, with contributions from SP, AT, PG, MFa, and LF. FA, SG, FD, and FS developed the S3M model, with evaluation contributions from SR, HS, EC, and UMdC. FA and SG validated IT-SNOW, with contributions from LR, SR, HS, EC, UMdC, AC, MFi, and OC. FP developed the historical dataset of the precipitation according to the modified conditional merging approach. GB processed streamflow data and assisted in the analysis of the streamflow–SWE relationship. LP, GS, and EF developed the operational chain providing daily Sentinel-2 snow-covered-area maps and assisted with snow-covered-area assimilation. FA carried out the analyses and prepared the paper, with input from all co-authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Financial support

This research has been supported by the Italian Civil Protection Department.

Review statement

This paper was edited by Baptiste Vandecrux and reviewed by three anonymous referees.


Alfieri, L., Avanzi, F., Delogu, F., Gabellani, S., Bruno, G., Campo, L., Libertino, A., Massari, C., Tarpanelli, A., Rains, D., Miralles, D. G., Quast, R., Vreugdenhil, M., Wu, H., and Brocca, L.: High-resolution satellite products improve hydrological modeling in northern Italy, Hydrol. Earth Syst. Sci., 26, 3921–3939,, 2022. a

Apicella, L., Puca, S., Lagasio, M., Meroni, A., Milelli, M., Vela, N., Garbero, V., Ferraris, L., and Parodi, A.: The predictive capacity of the high resolution weather research and forecasting model: a year-long verification over Italy, Bulletin of Atmospheric Science and Technology, 2, 1–14, 2021. a

Avanzi, F., Bianchi, A., Cina, A., De Michele, C., Maschio, P., Pagliari, D., Passoni, D., Pinto, L., Piras, M., and Rossi, L.: Centimetric Accuracy in Snow Depth Using Unmanned Aerial System Photogrammetry and a MultiStation, Remote Sens.-Basel, 10, 765,, 2018. a

Avanzi, F., Maurer, T., Glaser, S. D., Bales, R. C., and Conklin, M. H.: Information content of spatially distributed ground-based measurements for hydrologic-parameter calibration in mixed rain-snow mountain headwaters, J. Hydrol., 582, 124478,, 2020. a, b

Avanzi, F., Ercolani, G., Gabellani, S., Cremonese, E., Pogliotti, P., Filippa, G., Morra di Cella, U., Ratto, S., Stevenin, H., Cauduro, M., and Juglair, S.: Learning about precipitation lapse rates from snow course data improves water balance modeling, Hydrol. Earth Syst. Sci., 25, 2109–2131,, 2021a. a, b, c, d, e, f, g

Avanzi, F., Gabellani, S., Delogu, F., and Silvestro, F.: c-hydro/s3m-dev: (v5.1.0), Zenodo [software],, 2021b. a

Avanzi, F.: c-hydro/fp-s3m: (v1.0.1), Zenodo [software],, 2021c. a

Avanzi, F., Gabellani, S., Delogu, F., Silvestro, F., Cremonese, E., Morra di Cella, U., Ratto, S., and Stevenin, H.: Snow Multidata Mapping and Modeling (S3M) 5.1: a distributed cryospheric model with dry and wet snow, data assimilation, glacier mass balance, and debris-driven melt, Geosci. Model Dev., 15, 4853–4879,, 2022a. a, b, c, d, e

Avanzi, F., Gabellani, S., Delogu, F., Silvestro, F., Pignone, F., Bruno, G., Pulvirenti, L., Squicciarino, G., Fiori, E., Rossi, L., Puca, S., Toniazzo, A., Giordano, P., Falzacappa, M., Ratto, S., Stevenin, H., Cardillo, A., Fioletti, M., Cazzuli, O., Cremonese, E., Morra di Cella, U., and Ferraris, L.: IT-SNOW: a snow reanalysis for Italy blending modeling, in-situ data, and satellite observations, Zenodo [data set],, 2022b. a, b, c, d

Bales, R., Molotch, N. P., Painter, T. H., Dettinger, M. D., Rice, R., and Dozier, J.: Mountain hydrology of the western United States, Water Resour. Res., 42, W08432,, 2006. a, b

Barnett, T. P., Adam, J. C., and Lettenmaier, D. P.: Potential impacts of a warming climate on water availability in snow-dominated regions, Nature, 438, 303–309,, 2005. a

Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning Part I: numerical model, Cold Reg. Sci. Technol., 35, 123–145,, 2002. a

Bisci, C., Fazzini, M., Beltrando, G., Cardillo, A., and Romeo, V.: The February 2012 exceptional snowfall along the Adriatic side of Central Italy, Meteorol. Z., 21, 503–508,, 2012. a

Blöschl, G.: Scaling issues in snow hydrology, Hydrol. Process., 13, 2149–2175,<2149::AID-HYP847>3.0.CO;2-8, 1999. a

Bouamri, H., Boudhar, A., Gascoin, S., and Kinnard, C.: Performance of temperature and radiation index models for point-scale snow water equivalent (SWE) simulations in the Moroccan High Atlas Mountains, Hydrolog. Sci. J., 63, 1844–1862,, 2018. a

Bruno, G., Pignone, F., Silvestro, F., Gabellani, S., Schiavi, F., Rebora, N., Giordano, P., and Falzacappa, M.: Performing Hydrological Monitoring at a National Scale by Exploiting Rain-Gauge and Radar Networks: The Italian Case, Atmosphere, 12, 771,, 2021. a, b, c, d, e

Bruno, G., Avanzi, F., Gabellani, S., Ferraris, L., Cremonese, E., Galvagno, M., and Massari, C.: Disentangling the role of subsurface storage in the propagation of drought through the hydrological cycle, Adv. Water Resour., 169, 104305,, 2022. a

Bühler, Y., Adams, M. S., Bösch, R., and Stoffel, A.: Mapping snow depth in alpine terrain with unmanned aerial systems (UASs): potential and limitations, The Cryosphere, 10, 1075–1088,, 2016. a

Chiambretti, I., Dellavedova, P., Segor, V., Valt, M., and Cianfarra, P.: Winter 2013/2014 on the Italian Alps – Analysis and Lesson Learned About Avalanche Risk Treatment and Management Strategies, in: Proceedings of the ISSW 2014, International Snow Science Workshop 2014 Proceedings, Banff, Canada, (last access: 3 February 2023), 2014. a

Compo, G. P., Whitaker, J. S., Sardeshmukh, P. D., Matsui, N., Allan, R. J., Yin, X., Gleason, B. E., Vose, R. S., Rutledge, G., Bessemoulin, P., Brönnimann, S., Brunet, M., Crouthamel, R. I., Grant, A. N., Groisman, P. Y., Jones, P. D., Kruk, M. C., Kruger, A. C., Marshall, G. J., Maugeri, M., Mok, H. Y., Nordli, Ø., Ross, T. F., Trigo, R. M., Wang, X. L., Woodruff, S. D., and Worley, S. J.: The Twentieth Century Reanalysis Project, Q. J. Roy. Meteor. Soc., 137, 1–28,, 2011. a

Cox, L., Bartee, L., Crook, A., Farnes, P., and Smith, J.: The care and feeding of snow pillows, in: Proceedings of the 46th Annual Western Snow Conference, 18–20 April 1978, Otter Rock, Oregon, 40–47, 1978. a

Cui, G., Bales, R., Rice, R., Anderson, M., Avanzi, F., Hartsough, P., and Conklin, M.: Detecting Rain–Snow-Transition Elevations in Mountain Basins Using Wireless Sensor Networks, J. Hydrometeorol., 21, 2061–2081, 2020. a

Da Ronco, P., Avanzi, F., De Michele, C., Notarnicola, C., and Schaefli, B.: Comparing MODIS snow products Collection 5 with Collection 6 over Italian Central Apennines, Int. J. Remote Sens., 41, 4174–4205,, 2020. a, b

De Michele, C., Avanzi, F., Passoni, D., Barzaghi, R., Pinto, L., Dosso, P., Ghezzi, A., Gianatti, R., and Della Vedova, G.: Using a fixed-wing UAS to map snow depth distribution: an evaluation at peak accumulation, The Cryosphere, 10, 511–522,, 2016. a, b

Demirtaş, M.: The large-scale environment of the European 2012 high-impact cold wave: prolonged upstream and downstream atmospheric blocking, Weather, 72, 297–301,, 2017. a

DeWalle, D. R. and Rango, A.: Principles of Snow Hydrology, Cambridge University Press, Cambridge,, 2011. a

Dietz, A. J., Kuenzer, C., Gessner, U., and Dech, S.: Remote sensing of snow – a review of available methods, Int. J. Remote Sens., 33, 4094–4134,, 2012. a

Dozier, J.: Mountain hydrology, snow color, and the fourth paradigm, Eos T. Am. Geophys. Un., 92, 373–374,, 2011. a

Dozier, J., Bair, E. H., and Davis, R. E.: Estimating the spatial distribution of snow water equivalent in the world's mountains, WIREs Water, 3, 461–474,, 2016. a

Fiddes, J., Aalstad, K., and Westermann, S.: Hyper-resolution ensemble-based snow reanalysis in mountain regions using clustering, Hydrol. Earth Syst. Sci., 23, 4717–4736,, 2019. a, b

Fierz, C., Armstrong, R., Durand, Y., Etchevers, P., Greene, E., McClung, D., Nishimura, K., Satyawali, P., and Sokratov, S.: The International Classification for Seasonal Snow on the Ground, Tech. rep., IHP-VII Technical Documents in Hydrology N 83, IACS Contribution N 1, UNESCO – IHP, Paris, 2009. a

Flanner, M. G., Shell, K. M., Barlage, M., Perovich, D. K., and Tschudi, M. A.: Radiative forcing and albedo feedback from the Northern Hemisphere cryosphere between 1979 and 2008, Nat. Geosci., 4, 151–155,, 2011. a

Froidurot, S., Zin, I., Hingray, B., and Gautheron, A.: Sensitivity of Precipitation Phase over the Swiss Alps to Different Meteorological Variables, J. Hydrometeorol., 15, 685–696,, 2014. a

Gascoin, S., Grizonnet, M., Bouchet, M., Salgues, G., and Hagolle, O.: Theia Snow collection: high-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data, Earth Syst. Sci. Data, 11, 493–514,, 2019a. a

Gascoin, S., Grizonnet, M., Bouchet, M., Salgues, G., and Hagolle, O.: Theia Snow collection: high-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data, Earth Syst. Sci. Data, 11, 493–514,, 2019b. a

Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2), J. Climate, 30, 5419–5454,, 2017. a

Grünewald, T. and Lehning, M.: Are flat-field snow depth measurements representative? A comparison of selected index sites with areal snow depth measurements at the small catchment scale, Hydrol. Process., 29, 1717–1728,, 2015. a

Grünewald, T., Schirmer, M., Mott, R., and Lehning, M.: Spatial and temporal variability of snow depth and ablation rates in a small mountain catchment, The Cryosphere, 4, 215–225,, 2010. a

Harder, P., Schirmer, M., Pomeroy, J., and Helgason, W.: Accuracy of snow depth estimation in mountain and prairie environments by an unmanned aerial vehicle, The Cryosphere, 10, 2559–2571,, 2016. a

Harpold, A. A., Dettinger, M., and Rajagopal, S.: Defining snow drought and why it matters, EOS, 98,, 2017. a

Harrison, B. and Bales, R.: Skill Assessment of Water Supply Forecasts for Western Sierra Nevada Watersheds, J. Hydrol. Eng., 21, 04016002,, 2016. a, b, c

Hatchett, B. J. and McEvoy, D. J.: Exploring the Origins of Snow Drought in the Northern Sierra Nevada, California, Earth Interact., 22, 1–13,, 2018. a

Hatchett, B. J., Daudert, B., Garner, C. B., Oakley, N. S., Putnam, A. E., and White, A. B.: Winter Snow Level Rise in the Northern Sierra Nevada from 2008 to 2017, Water, 9, 899,, 2017. a

Huning, L. S. and AghaKouchak, A.: Global snow drought hot spots and characteristics, P. Natl. Acad. Sci. USA, 117, 19753–19759,, 2020. a

Immerzeel, W. W., van Beek, L. P. H., and Bierkens, M. F. P.: Climate Change Will Affect the Asian Water Towers, Science, 328, 1382–1385, 2010. a

Immerzeel, W. W., Lutz, A. F., Andrade, M., Bahl, A., Biemans, H., Bolch, T., Hyde, S., Brumby, S., Davies, B. J., Elmore, A. C., Emmer, A., Feng, M., Fernández, A., Haritashya, U., Kargel, J. S., Koppes, M., Kraaijenbrink, P. D. A., Kulkarni, A. V., Mayewski, P. A., Nepal, S., Pacheco, P., Painter, T. H., Pellicciotti, F., Rajaram, H., Rupper, S., Sinisalo, A., Shrestha, A. B., Viviroli, D., Wada, Y., Xiao, C., Yao, T., and Baillie, J. E. M.: Importance and vulnerability of the world's water towers, Nature, 577, 364–369, 2020. a

Kirchner, P. B., Bales, R. C., Molotch, N. P., Flanagan, J., and Guo, Q.: LiDAR measurement of seasonal snow accumulation along an elevation gradient in the southern Sierra Nevada, California, Hydrol. Earth Syst. Sci., 18, 4261–4275,, 2014. a

Kling, H., Fuchs, M., and Paulin, M.: Runoff conditions in the upper Danube basin under an ensemble of climate change scenarios, J. Hydrol., 424–425, 264–277,, 2012. a, b

Knoben, W. J. M., Freer, J. E., and Woods, R. A.: Technical note: Inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, Hydrol. Earth Syst. Sci., 23, 4323–4331,, 2019. a

Lagasio, M., Fagugli, G., Ferraris, L., Fiori, E., Gabellani, S., Masi, R., Mazzarella, V., Milelli, M., Parodi, A., Pignone, F., Puca, S., Pulvirenti, L., Silvestro, F., Squicciarino, G., and Parodi, A.: A Complete Meteo/Hydro/Hydraulic Chain Application to Support Early Warning and Monitoring Systems: The Apollo Medicane Use Case, Remote Sens.-Basel, 14, 14, 6348,, 2022. a

Li, D., Wrzesien, M. L., Durand, M., Adam, J., and Lettenmaier, D. P.: How much runoff originates as snow in the western United States, and how will that change in the future?, Geophys. Res. Lett., 44, 6163–6172,, 2017. a

Lievens, H., Demuzere, M., Marshall, H.-P., Reichle, R. H., Brucker, L., Brangers, I., de Rosnay, P., Dumont, M., Girotto, M., Immerzeel, W. W., Jonas, T., Kim, E. J., Koch, I., Marty, C., Saloranta, T., Schöber, J., and De Lannoy, G. J. M.: Snow depth variability in the Northern Hemisphere mountains observed from space, Nat. Commun., 10, 1–12, 2019. a, b, c, d, e

Liu, Y., Fang, Y., and Margulis, S. A.: High Mountain Asia UCLA Daily Snow Reanalysis, Version 1, NASA National Snow and Ice Data Center Distributed Active Archive Center,, 2021. a

Lundquist, J. D., Hughes, M., Henn, B., Gutmann, E. D., Livneh, B., Dozier, J., and Neiman, P.: High-Elevation Precipitation Patterns: Using Snow Measurements to Assess Daily Gridded Datasets across the Sierra Nevada, California, J. Hydrometeorol., 16, 1773–1792,, 2015. a

Magnusson, J., Wever, N., Essery, R., Helbig, N., Winstral, A., and Jonas, T.: Evaluating snow models with varying process representations for hydrological applications, Water Resour. Res., 51, 2707–2723, 2015. a

Main-Knorn, M., Pflug, B., Louis, J., Debaecker, V., Müller-Wilm, U., and Gascon, F.: Sen2Cor for Sentinel-2, in: Image and Signal Processing for Remote Sensing XXIII, edited by: Bruzzone, L., vol. 10427, p. 1042704, International Society for Optics and Photonics, SPIE,, 2017. a

Malek, S. A., Avanzi, F., Brun-Laguna, K., Maurer, T., Oroza, C. A., Hartsough, P. C., Watteyne, T., and Glaser, S. D.: Real-Time Alpine Measurement System Using Wireless Sensor Networks, Sensors, 17, 2583,, 2017. a

Margulis, S. A., Girotto, M., Cortés, G., and Durand, M.: A Particle Batch Smoother Approach to Snow Water Equivalent Estimation, J. Hydrometeorol., 16, 1752–1772,, 2015. a, b

Margulis, S. A., Cortés, G., Girotto, M., and Durand, M.: A Landsat-Era Sierra Nevada Snow Reanalysis (1985–2015), J. Hydrometeorol., 17, 1203–1221,, 2016. a, b

Montoya, E. L., Dozier, J., and Meiring, W.: Biases of April 1 snow water equivalent records in the Sierra Nevada and their associations with large-scale climate indices, Geophys. Ress Letts, 41, 1–7, 2014. a

Musselmann, K. N., Clark, M. P., Liu, C., Ikeda, K., and Rasmussen, R.: Slower snowmelt in a warmer world, Nat. Clim. Change, 7, 214–219,, 2017. a

Niwano, M., Suya, M., Nagaya, K., Yamaguchi, S., Matoba, S., Harada, I., and Ohkawara, N.: Estimation of Seasonal Snow Mass Balance all over Japan Using a High-Resolution Atmosphere-Snow Model Chain, SOLA, 18, 193–198,, 2022. a

Olefs, M., Koch, R., Schöner, W., and Marke, T.: Changes in Snow Depth, Snow Cover Duration, and Potential Snowmaking Conditions in Austria, 1961–2020—A Model Based Approach, Atmosphere, 11, 1330, 2020. a, b, c

Pagano, T., Garen, D., and Sorooshian, S.: Evaluation of Official Western U. S. Seasonal Water Supply Outlooks, 1922–2002, J. Hydrometeorol., 5, 896–909,<0896:EOOWUS>2.0.CO;2, 2004. a

Pagano, T. C., Wood, A. W., Ramos, M.-H., Cloke, H. L., Pappenberger, F., Clark, M. P., Cranston, M., Kavetski, D., Mathevet, T., Sorooshian, S., and Verkade, J. S.: Challenges of Operational River Forecasting, J. Hydrometeorol., 15, 1692–1707,, 2014. a

Painter, T. H., Berisford, D. F., Boardman, J. W., Bormann, K. J., Deems, J. S., Gehrke, F., Hedrick, A., Joyce, M., Laidlaw, R., Marks, D., Mattmann, C., McGurk, B., Ramirez, P., Richardson, M., Skiles, S. M., Seidel, F. C., and Winstral, A.: The Airborne Snow Observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and snow albedo, Remote Sens. Environ., 184, 139–152,, 2016. a, b

Pfeffer, W. T., Arendt, A. A., Bliss, A., Bolch, T., Cogley, J. G., Gardner, A. S., Hagen, J.-O., Hock, R., Kaser, G., Kienholz, C., Miles, E. S., Moholdt, G., Mölg, N., Paul, F., Radić, V., Rastner, P., Raup, B. H., Rich, J., Sharp, M. J., and The Randolph Consortium: The Randolph Glacier Inventory: a globally complete inventory of glaciers, J. Glaciol., 60, 537–552,, 2014. a

Piazzi, G., Thirel, G., Campo, L., and Gabellani, S.: A particle filter scheme for multivariate data assimilation into a point-scale snowpack model in an Alpine environment, The Cryosphere, 12, 2287–2306,, 2018. a

Pulliainen, J., Luojus, K., Derksen, C., Mudryk, L., Lemmetyinen, J., Salminen, M., Ikonen, J., Takala, M., Cohen, J., Smolander, T., and Norberg, J.: Patterns and trends of Northern Hemisphere snow mass from 1980 to 2018, Nature, 581, 294–298, 2020. a

Rice, R. and Bales, R. C.: Embedded-sensor network design for snow cover measurements around snow pillow and snow course sites in the Sierra Nevada of California, Water Resour. Res., 46, W03537,, 2010. a

Rolland, C.: Spatial and Seasonal Variations of Air Temperature Lapse Rates in Alpine Regions, J. Climate, 16, 1032–1046,<1032:SASVOA>2.0.CO;2, 2003. a

Rosenberg, E. A., Wood, A. W., and Steinemann, A. C.: Statistical applications of physically based hydrologic models to seasonal streamflow forecasts, Water Resour. Res., 47, W00H14,, 2011. a

Ryan, W. A., Doesken, N. J., and Fassnacht, S. R.: Evaluation of Ultrasonic Snow Depth Sensors for U. S. Snow Measurements, J. Atmos. Ocean. Tech., 25, 667–684,, 2008. a, b, c

Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., Behringer, D., Hou, Y.-T., Chuang, H.-Y., Iredell, M., Ek, M., Meng, J., Yang, R., Mendez, M. P., van den Dool, H., Zhang, Q., Wang, W., Chen, M., and Becker, E.: The NCEP Climate Forecast System Version 2, J. Climate, 27, 2185–2208,, 2014. a

Sinclair, S. and Pegram, G.: Combining radar and rain gauge rainfall estimates using conditional merging, Atmos. Sci. Lett., 6, 19–22, 2005. a

Slatyer, R. A., Umbers, K. D. L., and Arnold, P. A.: Ecological responses to variation in seasonal snow cover, Conserv. Biol., 36, e13727,, 2022. a

Soruco, A., Vincent, C., Rabatel, A., Francou, B., Thibert, E., Sicart, J. E., and Condom, T.: Contribution of glacier runoff to water resources of La Paz city, Bolivia (16 S), Ann. Glaciol., 56, 147–154,, 2015. a

Sturm, M. and Liston, G. E.: Revisiting the Global Seasonal Snow Classification: An Updated Dataset for Earth System Applications, J. Hydrometeorol., 22, 2917–2938,, 2021. a, b

Sturm, M., Holmgren, J., and Liston, G. E.: A seasonal snow cover classification system for local to global applications, J. Climate, 8, 1261–1283, 1995. a

Sturm, M., Goldstein, M. A., and Parr, C.: Water and life from snow: A trillion dollar science question, Water Resour. Res., 53, 3534–3544,, 2017. a, b

Tang, Q. and Lettenmaier, D. P.: Use of satellite snow-cover data for streamflow prediction in the Feather River Basin, California, Int. J. Remote Sens., 31, 3745–3762,, 2010. a

Terzago, S., Andreoli, V., Arduini, G., Balsamo, G., Campo, L., Cassardo, C., Cremonese, E., Dolia, D., Gabellani, S., von Hardenberg, J., Morra di Cella, U., Palazzi, E., Piazzi, G., Pogliotti, P., and Provenzale, A.: Sensitivity of snow models to the accuracy of meteorological forcings in mountain environments, Hydrol. Earth Syst. Sci., 24, 4061–4090,, 2020. a

Toreti, A., Bavera, D., Avanzi, F., Cammalleri, C., De Felice, M., De Jager, A., Di Ciollo, C., Gabellani, S., Maetens, W., Magni, D., Manfron, G., Masante, D., Mazzeschi, M., Mccormick, N., Naumann, G., Niemeyer, S., Rossi, L., Seguini, L., Spinoni, J., and Van Den Berg, M.: Drought in northern Italy – March 2022: GDO analytical report, European Commission, Joint Research Centre Luxembourg,, 2022.  a

Valery, A., Andréassian, V., and Perrin, C.: Inverting the hydrological cycle: when streamflow measurements help assess altitudinal precipitation gradients in mountain areas, IAHS Publ., 333, 281–286, 2009. a

Valt, M., Guyennon, N., Salerno, F., Petrangeli, A. B., Salvatori, R., Cianfarra, P., and Romano, E.: Predicting new snow density in the Italian Alps: A variability analysis based on 10 years of measurements, Hydrol. Process., 32, 3174–3187,, 2018. a

Vernay, M., Lafaysse, M., Monteiro, D., Hagenmuller, P., Nheili, R., Samacoïts, R., Verfaillie, D., and Morin, S.: The S2M meteorological and snow cover reanalysis over the French mountainous areas: description and evaluation (1958–2021), Earth Syst. Sci. Data, 14, 1707–1733,, 2022. a, b

Viviroli, D., Messerli, H. H. D. B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43, W07447,, 2007. a, b, c

WMO: Guide to Instruments and Methods of Observation, World Meteorological Organization (WMO), Geneva, 2018. a

Zanotti, F., Endrizzi, S., Bertoldi, G., and Rigon, R.: The GEOTOP snow module, Hydrol. Process., 18, 3667–3679,, 2004. a

Zhang, Z., Glaser, S., Bales, R., Conklin, M., Rice, R., and Marks, D.: Insights into mountain precipitation and snowpack from a basin-scale wireless-sensor network, Water Resour. Res., 53, 6626–6641,, 2017a. a, b

Zhang, Z., Glaser, S. D., Bales, R. C., Conklin, M., Rice, R., and Marks, D.: Technical Report: The Design and Evaluation of a Basin-scale Wireless Sensor Network for Mountain Hydrology, Water Resour. Res., 53, 4487–4498,, 2017b. a

Short summary
Snow cover has profound implications for worldwide water supply and security, but knowledge of its amount and distribution across the landscape is still elusive. We present IT-SNOW, a reanalysis comprising daily maps of snow amount and distribution across Italy for 11 snow seasons from September 2010 to August 2021. The reanalysis was validated using satellite images and snow measurements and will provide highly needed data to manage snow water resources in a warming climate.
Final-revised paper