Snow cover prediction in the Italian Central Apennines using weather forecast and snowpack numerical models

Italy is a territory characterized by complex topography with the Apennines mountain range crossing the entire peninsula with its highest peaks in central Italy. Using the latter as area of interest and the winter season during 2018-2019, the goal of this study is to investigate the ability of snow cover models to reproduce the observed snow height, using forecast weather data as meteorological forcing. We here consider two well-known ground surface and soil models: i) Noah LSM, a single-layer 5 Eulerian model; ii) Alpine3D, a multi-layer Lagrangian model. We adopt the Weather Research and Forecasting (WRF) model to produce the meteorological data to drive both Noah LSM and Alpine3D at regional scale with a spatial resolution of 3 km. While Noah LSM is already online coupled with the WRF model, we develop here a dedicated offline coupling between WRF and Alpine3D LSM. We validate the WRF simulations of surface meteorological variables in central Italy using a dense network of automatic weather stations, obtaining correlation coefficients of 0.84, 0.58, 0.4, 0.77 and 0.66 for air temperature, relative 10 humidity, wind speed, incoming shortwave radiation and daily precipitation, respectively. The performances of both WRFNoah and WRF-Alpine3D, are evaluated by comparing simulated and measured snow heights, provided by a quality-controlled network of snow stations located in Central Apennines. We find that WRF-Noah and WRF-Alpine3D models present similar correlation coefficients equal to 0.77 and 0.71, respectively, but the WRF-Alpine3D model produces a lower bias (about 2.2 cm) compared to the WRF-Noah model (about -8.0 cm) in the snow height estimation. For the estimation of daily snow height 15 variation WRF-Noah and WRF-Alpine3D present similar results with correlation coefficients of 0.72 and 0.71, respectively, but again WRF-Alpine3D showed a bias lower than WRF-Noah, about 0.09 cm and -0.22 cm respectively. The WRF-Noah model is slightly better than WRF-Alpine3D to reproduce the snow cover area observed with respect to the Moderate Resolution Imaging Spectroradiometer (MODIS) with the Jaccard spatial correlation index of 0.38 and 0.36 (optimal value equal 1), respectively, and Average Symmetric Surface Distance (ASSD) of 2.0 and 2.2 (optimal value equal 0), respectively, even 20 though both models tend to overestimate it. We finally show that snow settlement rate in WRF-Alpine3D is mainly driven by densification, whereas in WRF-Noah there is also an important contribution of snow melting especially at high elevation. As a general conclusion, the snow cover extension and height in central Italy at moderate spatial resolution (3 km) are well reproduced by both WRF-Noah and WRF-Alpine3D, but with the latter exhibiting a lower bias likely due to its multi-layer more sophisticated numerical scheme. 25 1 https://doi.org/10.5194/tc-2021-285 Preprint. Discussion started: 3 November 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Snow cover is a key element within global and local climate due to its high reflectivity of solar radiation and thermal insulation effects of land surfaces, inducing a non linear feedback in the Earth energy balance (Hall, 2004).Snowpack melting is an essential water supply to sub-surface reservoirs, urban populations and agricultural activities which mainly rely on frozen water stored in the snowpack (Barnett et al., 2005).For mountainous touristic places the persistence of snow cover is the basis of their winter sport economy (Vanat, 2020).Nevertheless, the snowpacl can be potentially dangerous when sliding down a slope causing avalanches with possible major damages to local flora, fauna and even human settlements (Bebi et al., 2009).
The simulation of the seasonal snow cover is particularly challenging in mountainous regions because of the complex interaction between atmospheric flow and topography.Indeed the synoptic flow approaching the topography may lead to the formation of local phenomena, like orographic precipitations, that have a large influence on snow deposition patterns and evolution (Roe, 2005).Also small scale phenomena, like wind erosion or accumulation (Sommer et al., 2018) and preferential deposition (Lehning et al., 2008), can cause large snow cover variations within few meters.Thus the quality of a distributed simulation of the snow cover is closely related to the horizontal resolution and the quality of the atmospheric data used to force the snow cover model.Recently the use of numerical weather prediction (NWP) models to drive snow cover models have become more and more investigated, thanks to the improved computer performances allowing to increase the spatial resolution and decrease the computational time.Bellaire et al. (2011), Bellaire et al. (2013) and Bellaire and Jamieson (2013) simulated the snow cover properties over British Columbia with the SNOWPACK model (Bartelt and Lehning, 2002;Lehning et al., 2002a, b), forced with atmospheric data provided by the the NWP model GEM15 (Mailhot et al., 2006), suggesting the possibility to force snow cover models with simulated weather data where observations are not available, and showing promising results on the possibility of predicting critical layers formation.Still in western Canada, Schirmer and Jamieson (2015) showed that forcing SNOWPACK with the NWP model GEM-LAM (Erfani et al., 2005) at 2.5 km resolution provided better results than using the GEM15 model, because of its higher resolution.Moreover, Horton et al. (2015) and Horton and Jamieson (2016) investigated the dependency of surface hoar formation on local meteorological and terrain conditions and the possibility to predict it driving SNOWPACK with the GEM-LAM model.Vionnet et al. (2016) simulated the snowpack evolution over the French Alps driving the snow cover model Crocus (Vionnet et al., 2012) with the NWP model AROME (Seity et al., 2011), and Quéno et al. (2016) applied the same model chain on French and Spanish Pyrenees.Bellaire et al. (2017) used SNOWPACK forced with the NWP model COSMO (Doms and Schättler, 2002) to forecast wet snow avalanche activity on the Swiss Alps.Vionnet et al. (2017) forced Crocus with the NWP model Meso-NH (Lafore et al., 1997) at 50 meter resolution, showing that for their case study in the French Alps, snow eolian transport is the main cause of snow height variability.Gerber et al. (2018) compared COSMO-WRF simulations (Skamarock et al., 2008) from 135 m to 50 m resolution over the Swiss Alps to radar estimations, and highlighted that in presence of complex terrain a good representation of the topography is essential to predict the observed snow precipitation and accumulation variability.Lastly, Sharma et al. (2021)  to WRF.This permits to simulate a large number of snow layer, and thanks to the blowing snow scheme even snow erosion and redistribution can be simulated, as they showed for case studies in Antarctica and Swiss Alps.
Most of the cited works are focused on higher latitude and altitude mountain ranges, like the Canadian mountains, European Alps or Pyrenees, which have many peaks above 3000 m, or on cold regions like Antarctica.So far no study focused on lower latitude mountain ranges with peaks below 3000 m, such as the Italian Apennines.The latter is a series of mountains, bordered by narrow coastal lands, forming a great arc along the Italian peninsula, from the Maritime Alps down to the Egadi Islands in western Sicily.Their total length is about 1400 km and their width goes from about 40 to 200 km.Their eastern slopes down to the Adriatic Sea are typically steeper and more verdant than the western ones.As a matter of fact, the Apennines play a significant role in the Italian climate being a natural topographic barrier to both western Atlantic cyclonic fronts and eastern cold air intrusions, causing intense snowfalls and hazardous avalanche activity (Chiambretti and Sofia, 2018), which may lead to fatal tragedies (Frigo et al., 2021).Its highest peak is the mount Corno Grande (2912 m) in the Gran Sasso chain in central Italy, belonging to the Central Apennines, in the Abruzzo region.The Gran Sasso is also hosting the southernmost European paraglacier (glacial apparatus), named Calderone, a key sentinel of current climate change in the Mediterranean region.Due to its complex topography and relatively low altitudes, the Apennines snow cover monitoring and forecast are quite cumbersome and, to some extent, open further issues with respect to high altitude mountainous regions.
The aim of this study is to investigate the ability of two snow cover numerical models to reproduce the observed snow height in Central Apennines, using forecast weather data as meteorological forcing.To this purpose, we use two well-known snow and soil models: i) Noah LSM, a single-layer Eulerian model (Chen and Dudhia, 2001); ii) Alpine3D, a multi-layer Lagrangian model (Lehning et al., 2006).We adopt the Weather Research and Forecasting (WRF) model to produce the meteorological data to drive both Noah LSM and Alpine3D at regional scale with a spatial resolution of 3 km.
While Noah LSM is already online coupled with the WRF model, we develop here a dedicated offline coupling between WRF and Alpine3D for the first time, at least to the best of the authors knowledge.The offline coupling of the two models give the opportunity of running just Alpine3D in the case WRF simulations have already been performed.This is particularly useful if different parametrizations have to be tested in Alpine3D only and not in WRF.Moreover, the offline coupling of WRF and Alpine3D may be easily used by the services that already uses WRF operationally, by just installing Alpine3D and the interfacing library, without the need of modifying the already present WRF configuration.Nevertheless, it is important to note that the offline coupling of Alpine3D and WRF can be extended to other NWP models, just making small modifications to the interfacing library.
By means of the two numerical models WRF-Noah and WRF-Alpine3D, we simulate the snowpack evolution over the Italian Central Apennines for the winter 2018-2019 and compare the snow cover forecast results with in situ measurements of snow height and satellite-based radiometric observations of snow extent.
The present paper is structured as follows.In Sect. 2 we provide a brief description of the climate of the area of interest and a synoptic view of winter 2018-2019.In Sect.3.1 we describe the observations data set, the atmospheric model and the two snow cover models, as well as the built numerical model chain.We then introduce the statistical indices used to evaluate the https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.model skills.Finally, in Sect. 4 we compare and discuss observed and simulated data in terms of: i) atmospheric forcing, ii) snow height, iii) snow cover extent.Conclusions are drawn in last section.

Area and period of interest
In this section we first introduce the climatology of Italian Central Apennines with a view on the synoptic characteristics of the winter 2018-2019.

Climatology of Italian Central Apennines
From a climatological point of view, central Italy is classified as Cfs, according to the Koppen system (Köppen, 1931;Belda et al., 2014).As a part of the Mediterranean region, the climate in this area of the Italian peninsula is characterized by alternating rainy periods in winter months and dry periods during Summer.In general, the thermal excursion, on daily and annual basis, does not exceed 21°C, due to the mitigating effect of both the Tyrrhenian and the Adriatic seas that surround a quite narrow strip of land, with an average width of 240 km in W-E direction.Besides, due to its complex topography and altitudes ranging from 0 to 2914 m a.s.l.(Corno Grande peak), in an extension of few tens of kilometres, the area is also characterized by a high micro-climatic variability, from meso-Mediterranean to sub-continental temperate, in the inner part, close to the highest peaks of the Apennines ridge (Lena et al., 2012).
According the statistics provided by Report no.55 of the Italian Institute for Environmental Protection and research (ISPRA, 2015) the annual average temperature in Central Italy ranged from 5°C (mountain areas) to 20°C (coastal area) in the reference period 1981-2010.The minimum and maximum temperature values are encompassed between 3°C and 11 °C and between 8°C and 22°C, respectively.Annual average accumulated precipitation is estimated to be between 600 and 1500 mm, being the yearly maxima mainly localized in the western slope of the Apennines ridge, exposed to the Tyrrhenian Sea.Minimum precipitation values occur over both the western and eastern littoral areas.Summer coincides with the dry period, when accumulated precipitation ranges between 40 to 80 mm.In this context, it is worth to mention as the average annual temperature and precipitation reported in the Central Apennines area are 3.7 °C (-8.9 °C the mean minimum temperature) and 1170 mm, respectively (Petriccione and Bricca, 2019).
In the same mountain area, the precipitation maxima are reached in Autumn and Spring and no drought period occurs during summer.As for the snow, mean solid precipitation depth (historical data from 1921 to 1960) on annual basis ranged from 0-5 cm on the Tyrrhenian littoral area to a maximum value 100-300 cm on the Apennines (Ministero dei lavori pubblici, 1973).
Eastern slope and Adriatic side was found to have a greater average amount of snow, with climatological values up to 20-50 cm in the coastal zone.Accordingly, the estimated number of chill days ranged from a minimum value of 0-1 in the western coast, to 150-200 days in the inner part (Ministero dei lavori pubblici, 1973).In the last decades, snowfall, snow cover and snow persistence over the area has decreased, although a few extreme snowfall events were recorded (Piacentini et al., 2020).

Synoptic view of winter 2018-2019
From a synoptic point view, winter 2018-2019 may be divided into three phases, as can be noticed from Fig. 1, where the 500 hPa geopotential height, 850 hPa temperature and precipitation monthly anomalies over Europe are displayed.During December, the Euro-Mediterranean area was characterized by a strong positive 500 hPa geopotential height anomaly centered on the Iberian Peninsula, whilst a near neutral, or slightly negative, geopotential anomaly affected Eastern Europe.This resulted in a warm and dry pattern over the Western Mediterranean basin.By contrast, in January an oceanic anticyclonic system with above normal geopotential heights determined a flux of Arctic cold air toward the Mediterranean through Eastern Europe.This generated a generalized cold temperature anomaly over Central and Eastern Europe and Central Mediterranean basin, with associated above normal precipitations over Central and Eastern Mediterranean Sea.Finally, February was dominated by a robust 500 hPa geopotential height anomaly extended to the whole continental Europe.As a result, the Euro-Mediterranean region was characterized by warm temperatures and dry conditions.This pattern affected the precipitations in central Italy as follows.In the second decade of December, the circulation over Central Italy was affected by geopotential negative anomalies centred on Eastern Europe.In the middle of December a series of small cold Northern or Northern-Eastern impulses leaded to snowfall down to 700 m a.s.l..During the second part of the month, the atmospheric conditions were influenced by the positive geopotential anomaly localized on the Iberian Peninsula.
The associated mild weather caused a consistent retreat of the snow cover.The most important snowstorm of the season occurred between 2 nd and 4 th January 2020, when an elevation of the Azores anticyclone over the British Isles determined a flux of Arctic maritime air on the Central Mediterranean Sea (see Fig. 2a and Fig. 2d).Minor snowfall events occurred till 23rd January, when a new impulse of cold air from North Atlantic, leaded to the formation of a low pressure system on the Mediterranean Sea (see Fig. 2b and Fig. 2e), resulting in a new snowfall on central Italy.Between 1st and 3rd February, the passage of an oceanic through drawn air from North Africa (see Fig. 2c and Fig. 2f), advecting warm air over Italy.The rest of the month was characterized by mild and dry atmospheric conditions.

Data, models and evaluation methods
In this section we present the observational data, we describe the considered atmospheric and snow cover models and finally we define the performances evaluation methods.

Observational data set
The observational data set consists of in situ measurements of atmospheric and snow cover conditions, coming from automatic weather stations (AWS) as well as of snow cover extent derived from satellite-based observations.The in situ measurements come from a dense network of 702 automatic weather stations (AWS), shown in Fig. 3b shortwave radiation, precipitation and snow height.Table 1 shows the number of AWS providing the above variables.The AWS have an acquisition interval ranging between 15 minutes and 1 hour: all observed variables are daily averaged, with the exception of precipitation which is accumulated over one day, in order to set a common temporal framework.The snow height measurements are carried out with automatic ultra-sonic sensors, installed on 13 weather stations and located along the Italian Central Apennines (see Fig. 3b).Identification number (ID), coordinates and altitude of each automatic snow station are shown in Table 2.
We used 90 maps of snow cover fraction, from 1 December 2018 to 28 February 2019, at 500 m resolution.We first process MODIS maps in order to obtain cloud-free snow cover fraction maps as follows: the pixels classified as clouds are masked, and for each cloud pixel a linear interpolation between previous and consecutive images is carried out where the pixels are not cloud covered, in order to infer their fractional snow cover area.Then we bilinearly interpolate the cloud-free snow cover fraction maps on a regular grid of resolution 3 km.
In this section, we first introduce the weather forecast model, adopted and tuned to the study area of this work, as well as the two snow cover models, specifically Noah LSM and Alpine3D.In the methodological paragraph we describe the error indices used to quantitatively evaluate the error statistics when comparing both WRF-Noah and WRF-Alpine3D with observational data.

Description of numerical models
The atmospheric and snow cover numerical models are here briefly described, highlighting the space-time set up and main physical parameterizations employed in this study.

Weather Research and Forecasting (WRF) model
The Weather Research and Forecasting (WRF) model is used in this work to describe the atmospheric evolution at mesoscale.WRF is a regional non-hydrostatic meteorological model including several options for the parameterizations of atmospheric physical processes (Skamarock et al., 2008).As shown in Fig. 3a, the WRF atmospheric model is configured with three twoway nested domains.The largest domain has a 27 km spatial resolution and the second domain has 9 km spatial resolution, whereas the inner domain has a cloud-resolving module at 3 km resolution.The vertical grid is constituted by 33 vertical levels extending from surface up to 50 hPa (about 20 km above the sea level) with the first level at 10 m from the surface.
The WRF main parameterizations include: i) the Yonsei University (YSU) scheme for the planetary boundary layer (Hong et al., 2006); ii) the Rapid Radiative Transfer Model for General circulation model (RRTMG) for visible and infrared radiation (Iacono et al., 2008); iii) the Thompson scheme for cloud microphysics (Thompson et al., 2008); iv) the Grell-Freitas for cumulus convection (Grell and Freitas, 2014) with exception for the cloud-resolving domain, where no cumulus parameterization is activated.Initial conditions for WRF were taken from 6-hourly operational NCEP analyses at a resolution of 27 km, whilst for the nested domains, boundary conditions are derived from domain 1 and 2 simulations, respectively.In order to improve the meteorological simulations, a spectral nudging toward NCEP analysis of wind, temperature, and geopotential is applied in https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License. the outer domain above 850 hPa.With exception of the first simulation of this series, the land surface and soil properties are restarted by using the previous meteorological prediction in order to avoid a discontinuous simulation of the snowpack in the study area.We have chosen this setup for WRF, because it resulted the best one to simulate the main observed meteorological variables according to many preliminary sensitivity tests that we carried out.

Noah Land Surface Model
Within the WRF model, land surface processes are simulated with the Noah LSM (Chen and Dudhia, 2001).Noah LSM is model that simulates skin temperature, snow height, snow water equivalent, snow density, soil moisture, soil temperature and surface energy balance.The snowpack is represented as a single layer of varying thickness, instead the soil is discretized in four layers of thickness 10, 30, 60 and 100 cm respectively.Noah LSM is online coupled with WRF, and it provides to the parent model surface and latent heat fluxes, and upward shortwave and longwave radiation.The improvements in the Noah snowpack model introduced by Livneh et al. (2010) and Wang et al. (2010) were taken into account in our simulations.
The WRF-Noah model is run from 1 December 2018 till 28 February 2019, since after the end of February 2019 the observed snow cover in the area of interest was almost negligible.For these three months, a series of 60-h simulations are performed starting at 12 UTC, with the first 12 h discarded as model spin up.

Alpine3D snow and soil model
Alpine3D is a three-dimensional snow cover, land surface and soil numerical model (Lehning et al., 2006).It includes the one-dimensional model SNOWPACK (Bartelt and Lehning, 2002;Lehning et al., 2002a, b) to simulate snow, land surface and soil properties, the SnowDrift module to simulate the wind transport of snow, the EBalance module to compute radiation fields taking into account topographic shading and reflections effects and a runoff module.SNOWPACK is the core of Alpine3D, it is developed using a Lagrangian scheme which permits to provide detailed description of snow stratigraphy, mass and energy balance, simulating up to 50 snow layers.In order to run an Alpine3D simulation, it is necessary to provide the initial state of snow and soil (if soil simulation is enabled), a digital elevation model, a land cover model and a meteorological input.
Alpine3D is not coupled with WRF so that it is necessary to interface the two models.The digital elevation model and the land cover model, used in the Alpine3D simulation, are the same extracted from the WRF model.We have enabled in Alpine3D the soil simulation, creating 4 soil layers of 10, 30, 60 and 100 cm respectively, and initialized them using the soil condition extracted from WRF.Since on the beginning of December 2018 there was no snow cover in the study domain, Alpine3D is initialized with snow-free conditions on each grid cell.To drive Alpine3D, we use the air temperature, relative humidity, wind speed, incoming shortwave radiation, incoming longwave radiation, total precipitation, precipitation phase and ground surface temperature extracted from WRF. WRF-Alpine3D hourly maps of snow height and snow water equivalent are then generated for the whole period.A scheme of the described WRF-Alpine3D model chain is shown in Fig. 4.
In our Alpine3D simulation setup, we turn off the SnowDrift and the EBalance modules: indeed, at the resolution of 3 km the wind transport of snow and the topographic shading and reflections effects have a negligible impact on the simulation results.
The canopy representation is also disabled.After a sensitivity test (not shown here), neutral atmospheric stability conditions https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.have been assumed as well as the selection of the "Zwart" and "Lehning new" parametrizations for new snow density and snow albedo respectively.Even if more computational demanding, we use the Richards water transport scheme instead of the more simple bucket scheme, because it improves the quality of the simulation results (Wever et al., 2014).

Performance evaluation indices
We here distinguish the error indices for atmospheric variable and snow height prediction from those used to estimate the discrepancy in the snow cover extension.
WRF skill in predicting the meteorological variables is evaluated comparing the output with in situ observations.The simulated air temperature, relative humidity, wind speed, incoming shortwave radiation and precipitation are bilinearly interpolated on the coordinates of the observed data.They are daily averaged, with the exception of precipitation rate which is accumulated over one day.The simulated variables are statistically evaluated using mean bias error (MBE), mean absolute error (MAE) and Pearson correlation coefficient (R).These are defined as follows: where F i and O i are the forecasted and observed values for each i of N observations and F and Ō the corresponding means.
The simulated snow heights are bilinearly interpolated on the coordinates of the in situ observations, and then re-sampled at a temporal resolution of one day.Starting from daily snow height, we derive the corresponding daily variations.In order to evaluate the model-based simulations of snow height and snow-height variation, we use the coefficients defined in Eq. ( 1), Eq.
(2), Eq. ( 3).For the snow height variation, the Equitable Threat Score index (ETS) is also used (as proposed by Nurmi (2003) and Schirmer and Jamieson (2015)), and quantify the model skill in hitting a binary observation, like snow height exceeding a certain threshold.The ETS is defined as: where HI (hits), FA (false alarms) and MI (misses) are defined in the contingency Table 3 and HI rdm is defined as follow: where α s = 2.6 is a distribution shape parameter, W s is the simulated snow water equivalent, W max is the snow water equivalent threshold above which the soil is 100% covered with snow, and varies locally according to MODIS land use classification.
By applying a threshold of 51%, we can build binary maps which indicate: i) snow absence if snow cover fraction is below the threshold; ii) snow presence if snow cover fraction area was equal or above the threshold.These binary maps are compared by using the Jaccard index (J) and the average symmetric surface distance (ASSD).J is defined as: where A and B are the binary matrices of the simulated and retrieved snow cover area maps.If the matrices perfectly overlap, J is equal to 1, instead if the matrices don't overlap J is equal to 0. To define the ASSD, we need first to introduce the modified directed Hausdorff distance (Dubuisson and Jain, 1994): where L A and L B are the boundaries of the matrices A and B, N L A is the number of elements in the boundary L A and d(a, L B ) is the Euclidean distance between point a and the closest point of boundary L B : The ASSD can be expressed as: and it ranges from 0 to ∞, where 0 means perfect overlap of the boundaries L A and L B .From the snow cover binary maps of MODIS, WRF-Noah and WRF-Alpine3D we also calculate the respective snow cover area fraction, which is defined as the number of cells classified with snow divided by the total number of cells of the study domain.We compare the simulated snow cover area fractions with the one obtained from MODIS and we evaluate them in terms of MBE, MAE and R. We also build

Results and discussion
Both models, WRF-Noah and WRF-Alpine3D, are compared with observational data for the study area for winter 2018-2019.
The performance analysis is carried out in terms of the error indices, defined in the previous section, for both atmospheric and snow cover forecast.

Performance analysis of the atmospheric model
As a first validation, Figure 5 shows the scatter plots of the comparison between simulated and measured daily air temperature, relative humidity, wind speed and incoming shortwave radiation and daily precipitation.The colors represent the altitude of the sites where the measurements were taken.The statistical scores, obtained from this intercomparison, are reported in Table 4 for the three elevation ranges defined in Section 3.3 and for all altitudes.
WRF presents for the air temperature at all elevation ranges a correlation coefficient higher than 0.8, with an overall score of 0.84 (Fig. 5a).The WRF simulations exhibit a slightly negative MBE, especially for elevations between 800 m and 1600 m, where it also shows the highest MAE.
The relative humidity presents a correlation coefficient higher than 0.7 for mid and high elevation bands but lower (0.54) for low elevations, with and overall score of 0.58 (Fig. 5b).The MBE is negative for low elevations and increases with the altitude.
The MAE is similar for low and mid elevations and increases at high elevations.
From Fig. 5c we note that that WRF drastically underestimates the wind speed at high elevations.Indeed at low elevations WRF shows a positive MBE, which decreases and becomes negative as the elevation increases.Clearly, also the MAE show an elevation dependence, and increases from low to higher altitudes.Moreover, the correlation coefficient is minimum at high elevations and maximum at mid elevations.The wind speed underestimation at high elevation is likely due to an underestimated elevation of the topography.Indeed the 3 km resolution of the model simulation smooths the highest peaks, leading to a lower simulated wind speed at station location.
Figure 5d shows that WRF overestimates the incoming shortwave radiation.At high elevations, WRF has the highest MAE and the lowest correlation with the in situ observations.The general overestimation of incoming shortwave radiation may be again partly due to the limited model.Indeed shading effect that influence the measured solar radiation may not be captured by the model, in which the topography is more smooth and low compared to reality. Figure 5e highlights that WRF reproduces the daily precipitation with a correlation coefficient between 0.67 (for low elevations) and 0.63 (for high elevations).WRF shows the tendency to overestimate the precipitation as the altitude increases.
Indeed, the bias is negligible below 800 m of altitude, but MBE reaches about 4 mm at the high altitude sites.This overestimation could be due to the fact that the rain gauge typically loses a fraction precipitation.When the precipitation is solid the underestimation is even larger because part of it is lost with sublimation, especially when the rain gauge is not heated.
In summary, for the WRF-simulated variables the overall correlation coefficient is above 0.5, except for the wind speed, because as we already discussed, the strong negative bias at high elevations decreases the overall correlation to 0.4.On the other hand, it is well known that the wind speed field is particularly variable in complex orographic regions and a point-like in situ measurement cannot be representative of 3 km spatially averaged wind speed at WRF model scales considered in this study.

Snow height
The time series intercomparison can be useful to show the temporal behaviour of the two models with respect to in situ observations.Figure 6a shows three the time series of daily means of observed and simulated snow height averaged over all available stations, reported in Table 1.Observations are indicated by a dashed gray line, while the numerical simulations are indicated by solid blue and red lines, corresponding to WRF-Noah and WRF-Alpine3D respectively.The shaded areas correspond to an interval of ±1 standard deviation.This comparison indicates that the snow height, simulated by the two models, is similar until mid January 2019, but during the last part of the simulation period, WRF-Alpine3D slightly overestimates the observed snow height, while WRF-Noah underestimates it.From mid January until the end of the simulation period, WRF-Noah also shows a smaller standard deviation when compared to WRF-Alpine3D, which means that the simulated values underestimate the observations in most of the 13th sites, and that WRF-Alpine3D has a better skill to capture the observed snow height variability.On the last month and an half of simulation the settlement rate simulated with WRF-Alpine3D is more similar to the observations compared to WRF-Noah.Another comparison of simulated and observed data is shown in the scatter plot of Fig. 7a, with blue and red dots indicating WRF-Noah and WRF-Alpine3D respectively.The plot shows that the regression lines for WRF-Noah and WRF-Alpine3D have almost the same slope (0.50 and 0.54 respectively), but the line corresponding to WRF-Noah has a lower intercept compared to WRF-Alpine3D (8.1 and 17.0 respectively).The correlation with the observations is slightly higher for Noah-LSM (0.77) compared to Alpine3D (0.71), however WRF-Noah has a negative bias compared to WRF-Alpine3D (-8.0 and 2.2 respectively), as can be seen from Table 5, which summarises the statistical indices of the comparison.
Snow height variation is an important indicator of the models skill to catch the dynamical behavior of the snow cover.
Figure 6b shows the time series of observed and simulated snow height variation, averaged over all available stations reported in Table 1.Observations are indicated by a solid gray line, while the simulations are indicated by solid blue and red lines, The plot shows that the regression lines for WRF-Noah and WRF-Alpine3D are almost equal, with similar slope (0.59 and 0.57 respectively) and intercept (-0.1 and 0.2 respectively).Also the correlation coefficient is similar between the two models (0.72 for WRF-Noah and 0.71 for WRF-Alpine3D), however some differences can be seen in the bias (Table 5), which is negative for WRF-Noah (-0.22) and positive for WRF-Alpine3D (0.09).
A more detailed view of the skill of WRF-Noah and WRF-Alpine3D in calculating the observed daily snow height variation comes from the analysis of Fig. 8a and 8b. Figure 8a shows the observed and modelled frequency distribution of daily snow height variation at the 13 validation sites used in this study.Grey column represents the observed variations and the blue and red columns represent the snow height variations simulated with WRF-Noah and WRF-Alpine3D, respectively.For the range of variations between -10 cm and +20 cm, WRF-Noah and WRF-Alpine3D show similar performances, simulating a number of snow height variations for each class similar to the observations.Outside this range, the behaviour of the two models is clearly different.Below -10 cm, WRF-Alpine3D shows a better skill in reproducing the observed snow height variation, especially in the interval [-30 cm, -20 cm], which is totally missed from WRF-Noah.WRF-Alpine3D performs better than WRF-Noah also in the range [20 cm, 30 cm], but its skills rapidly decrease for snow height variations larger than 30 cm, where WRF-Noah seems to be more accurate to reproduce the observations.The snow height variations observed in the range [40 cm, 50 cm] are not captured by both WRF-Noah and WRF-Alpine3D.Therefore, both models can well predict the central tendency of the distribution of the daily snow height variation, whereas they have more difficulties on the tails of the distribution.A further confirmation of this behaviour comes from Fig. 8b, where the ETS index for different snow thresholds is reported.A look at Fig. 8b shows again that WRF-Noah and WRF-Alpine3D has similar scores between -5 cm ad 20 cm thresholds, but they perform differently outside of this range.Indeed, compared to WRF-Alpine3D, WRF-Noah shows higher performances for snow height variations larger than 30 cm, however both models do not capture snow height variations larger than 40 cm.Instead, for variations smaller than -5 cm we can see that WRF-Alpine3D reaches scores larger than WRF-Noah.This highlights that WRF-Alpine3D has slightly better skills compared to WRF-Noah in predicting the observed snow settlement, and it is likely due to the more complex multi-layer Lagrangian scheme of the Alpine3D model, which guarantees a representation of the snowpack properties more sophisticated than the single-layer Eulerian Noah model, which is known to anticipate the complete melt of the snowpack even at 3 km resolution, as Pavelsky et al. (2011) show for Sierra Nevada Mountains, California.

Snow cover extent
MODIS satellite snow product is an essential independent data set to evaluate the skill of WRF-Noah and WRF-Alpine3D to reproduce the observed snow cover area.As already described, MODIS sensor is providing the snow cover fraction (SCF) maps with 500 m resolution.Thanks to the applied cloud-removal algorithm, we have been able to compare the snow cover extent also for the days with high cloud cover fraction.The snow cover area fraction (SCAF) can be computed over the chosen domain for each day of the study period, dividing the area classified as covered with snow by the total land area.The result of this calculation is reported in Fig. 12 and statistical indices obtained from this comparison are summarized in Table 6.Both models reproduce SCAF with a high correlation coefficient (0.93 and 0.90 for WRF-Noah and WRF-Alpine3D, respectively), but they overestimate the observed SCAF of more than 0.06%.WRF-Noah exhibits a slightly better performance with respect to WRF-Alpine3D in the estimation of the SCAF.
Finally, for each day of the study period, we have also evaluated J and ASSD spatial indices, comparing the SCA maps derived for WRF-Noah and WRF-Alpine3D with MODIS the ones derived from MODIS.J and ASSD time series are reported in Fig. 13 and their average values are summarized in Table 7.For J and ASSD WRF-Noah exhibits mean values of 0.38 and 2.0 respectively, whereas WRF-Alpine3D 0.36 and 2.2.This means that WRF-Noah has a slightly higher skill in reproducing the snow cover extent and shape observed with MODIS.

Sensitivity analysis of snow cover models
A sensitivity analysis of WRF-Noah and WRF-Alpine3D can highlight the main differences between the two snow models and give insights to improve the snow cover forecast.
Figure 14 shows WRF-Noah and WRF-Alpine3D maps of mean and maximum snow height, snow water equivalent, and the differences between the two models for each of these variables for the chosen study period.The comparison between WRF-Noah and WRF-Alpine3D mean snow height reveals that WRF-Alpine3D simulates a snowpack in average thinner compared to WRF-Noah on the upper part of the Apennines range, and thicker on the lower part.The difference between the two models in the upper part of the domain is likely due to shorter snow cover duration in WRF-Alpine3D compared to WRF-Noah, as it can be observed from Fig. 11d, which decreases the WRF-Alpine3D mean snow height over the simulation period.By contrast, in the lower part of our domain, the difference is not due to a longer snow cover duration in WRF-Alpine3D compared to WRF-Noah for that areas (Fig. 11d), and it is likely caused by a faster snowpack settlement in WRF-Noah which results in a thinner snowpack.The maps of maximum snow height are consistent with the maps of the mean snow height, except for the southern-east part of the domain, where WRF-Noah simulates a thicker snowpack when compared to the WRF-Alpine3D one.
This can be likely due to some snowfall events which brought less dense snow in WRF-Noah with respect to WRF-Alpine3D in that zone, causing a thicker snowpack, for just one or few days, but not in average for the entire study period.
The comparison of the maps of the mean and maximum snow water equivalent (SWE) reveal other important differences between the two models.The map of mean snow water equivalent (Fig. 14k) suggests that more mass is stored in average in the WRF-Noah snowpack with respect to WRF-Alpine3D for the upper part of the Apennines range, while less mass is stored in the lower part.As discussed about the mean snow height, the reason of the differences between the two models in the upper part of the domain is probably caused by a shorter snow cover duration in WRF-Alpine3D compared to WRF-Noah.
The differences in the lower part, instead, are not due to a longer snow cover duration in WRF-Alpine3D with respect to WRF-Noah, as we already discussed, and it may be not caused by larger snowfalls simulated by WRF-Alpine3D compared to WRF-Noah, because the input precipitation and temperature fields, and thus the solid precipitations, are the same in both models.The cause is likely due to a faster snow melt in WRF-Noah compared to WRF-Alpine3D, which reduces the mass stored in the snowpack through melt water percolation and runoff.The map of maximum snow water equivalent (Fig. 14l) is consistent with the map of mean snow water equivalent and it is worth nothing that in the southern-east part of the domain is not present the blue area visible in Fig. 14j.This note supports the hypothesis of a sporadic event of less dense fresh fallen snow in WRF-Noah compared to WRF-Alpine3D, which causes for WRF-Noah a peak in the snow height in that zone, but not in the snow water equivalent.
In order to investigate more in details the differences between WRF-Noah and WRF-Alpine3D in the temporal evolution of the snowpack at the elevation range intervals defined in Section 3.3, we have calculated for each elevation range the daily spatial averages of snow height and snow water equivalent.Figure 15 shows the resulting time series for both WRF-Noah (dashed line) and WRF-Alpine3D (solid line) at different elevation bands.At high elevations, the snow height accumulation was larger for WRF-Noah than WRF-Alpine3D, whilst both models simulated similar values for snow water equivalent accumulation.This behavior can be caused by a less dense fresh fallen snow in WRF-Noah compared to WRF-Alpine3D, which results in a thicker new snow layer, or on the other hand, by a sudden compaction of the fresh fallen snow due to overburden stress in the snowpack simulated by WRF-Alpine3D, that in turn results in a thinner snowpack compared to WRF-Noah.This difference is less evident between 800 m and 1600 m, and it seems to be absent at lower elevations.Is is also evident that for mid and high elevation the WRF-Alpine3D simulates a much more persistent snow cover compared to WRF-Noah, especially for the elevations above 1600 m, where the snow water equivalent in WRF-Alpine3D increases almost monotonically, while it drastically decreases in WRF-Noah in the last part of the study period.This means that in WRF-Noah simulation, the snowpack at high altitude is subject to melting, water percolation and subsequent runoff, which causes a loss of mass of the snowpack.This implies also a faster reduction of the snow height in WRF-Noah compared in WRF-Alpine3D, suggesting that in WRF-Alpine3D the main driver of the snowpack settlement is the snow densification, while in WRF-Noah there is also a large contribution of the snow melt.

Conclusions
In this paper we have shown the different skills of Noah LSM and Alpine3D models to simulate the snow cover in Italian Central Apennines when forced with WRF atmospheric data during winter 2018-2019.The study area is novel to snow cover studies, since most of past works focused on higher latitude and altitude mountain ranges or on cold regions.So far no study focused on lower latitude mountain ranges with peaks below 3000 m, like Central Apennines.
The performances of the WRF model to simulate air temperature, relative humidity, wind speed, incoming shortwave radiation and precipitation have been first evaluated, comparing them observations derived from a dense network of automatic weather stations.The WRF prediction model is capable to predict the observed atmospheric variables with a correlation coefficient always higher than 0.58, except for the wind speed for which the model bias increases with the altitude, likely because of the underestimation of the model topography, caused by the limited spatial resolution of the simulation.We have also showed that the snow height, simulated by WRF-Noah, has a bias of about -8.0 cm, whilst WRF-Alpine3D presents a bias of about 2.2 cm, thus smaller in the absolute value.In the daily snow height variation the two models exhibits a bias of opposite sign, WRF-Noah bias being of -0.22 cm and WRF-Alpine3D bias of 0.09 cm (thus also in this case smaller than WRF-Noah in the absolute value).WRF-Noah and WRF-Alpine3D have similar skills to predict small daily snow height variations, but the former has slightly better performance to reproduce daily variation larger than 30 cm, while the latter has higher skills to predict variation smaller than -10 cm.Indeed, the sophisticated Largangian scheme of Alpine3D permits to simulate snow settlement more accurately, especially in melting conditions.
Both models WRF-Noah and WRF-Alpine3D tend to overestimate the snow cover area fraction, provided by satellite-based MODIS products.WRF-Noah reveals to perform slightly better than WRF-Alpine3D with a bias of 0.063% as compared to 0.066%.The evaluation of both models performance with the Jaccard index (J) and the Average Symmetric Surface Distance (ASSD) confirmed the slightly better abilities of WRF-Noah to predict shape, location and extension of the observed snow cover.

List of FIGURES
online coupled WRF and SNOWPACK, and introduced a new blowing snow scheme.In this configuration, WRF drives SNOWPACK which acts as land surface model and gives feedback https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.
, maintained by the Italian Civil Protection (Italian Civil Protection Department and CIMA Research Foundation, 2014) and covering the entire study area domain.The measured variables are surface air temperature, relative humidity, wind speed, incoming https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.
5) https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.with N = HI + FA + MI + CR with CR standing for correct rejection.The skill of WRF-Noah and WRF-Alpine3D model in reproducing the snow cover extent observed with MODIS are also evaluated.To this end, snow cover fraction maps are obtained for both WRF-Noah and WRF-Alpine3D from the following equation, equivalent to the empirical relation found by Koren et al. (1999): https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.andcompare snow cover duration maps for MODIS, WRF-Noah and WRF-Alpine3D counting the number of times that a cell was classified as covered with snow in the binary maps.Then we compare WRF-Noah and WRF-Alpine3D maps of mean and max seasonal snow height and snow water equivalent in order to further investigate the causes of the differences between the two models.Finally, to further discern how the two models behave at different elevations, the study domain is classified into three elevation ranges (low elevation: <800 m, mid elevation: 800-1600 m, high elevation: ≥1600 m), we calculate for each range the daily spatial averages of snow height and snow water equivalent, obtaining the corresponding time series.
corresponding to WRF-Noah and WRF-Alpine3D respectively.The 2 models have similar skill in reproducing the observed https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.snow height variation during December and January, where most of the accumulation events occurred.During February, when only few episodes of solid precipitation occurred and the snow settlement accelerated, WRF-Alpine3D shows a higher skill compared to WRF-Noah in predicting the observed snow height variation.We further compared the simulated and observed daily snow height variation as shown in Fig. 7b, with blue and red dots indicating WRF-Noah and WRF-Alpine3D respectively.
Figure10shows the forested area obtained from the aggregation of the broad leaved forest, coniferous forest, and mixed forest classes of CORINE land use classification at 100 m horizontal resolution, reprojected at 3 km horizontal resolution.It cannot be excluded that WRF-Noah and WRF-Alpine3D overestimation of the snow cover area of Fig.9is partially due to a MODIS underestimation of the snow cover area where forests are present, as highlighted byGascoin et al. (2015).Duration of the snow cover is another important parameter, which has direct implications on many aspects, like local climate, water supply and flora and fauna development.Figure11shows the observed and modelled snow cover duration maps, obtained counting the number of times that each cell is classified as covered with snow within MODIS, WRF-Noah and WRF-Alpine3D derived SCA maps.Figure11also shows the differences between the snow cover duration simulated with WRF-Noah and WRF-Alpine3D and the one obtained from MODIS.The difference between WRF-Alpine3D and WRF-Noah snow duration is also reported.These figures suggest that WRF-Noah and WRF-Alpine3D tend to simulate a longer snow cover duration on the entire Central Apennines mountain range.The best agreement of the numerical models with the satellite-based observations is found on the highest peaks and at the valley bottoms.By contrast, the largest differences emerge within middle elevation zones, where both WRF-Noah and WRF-Alpine3D simulate a much longer snow cover duration.A closer inspection of figure 11d highlights noticeable differences between WRF-Noah and WRF-Alpine3D in the snow cover persistence.WRF-Noah simulates a more persistent snow cover when compared to WRF-Alpine3D over most of the Apennines mountain range, with the exception over the Abruzzo region in the central-east part of the domain, where WRF-Alpine3D predicts a longer snow cover duration.
https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.Finally, we have intercompared WRF-Noah and WRF-Alpine3D in terms of snow height and snow water equivalent, showing that, at high elevations, WRF-Alpine3D simulates a thicker mean seasonal snowpack compared to WRF-Noah.This highlights that for WRF-Alpine3D the snow settlement is mainly driven by compaction, whilst in WRF-Noah an important contribution of melting can be also observed, which accelerates the snow settlement.Future work should be oriented to gather further in situ measurements of snow water equivalent, density, temperature, grain shape and grain size, necessary to quantify the abilities of WRF-Noah and WRF-Alpine3D to reproduce the observed snow settlement.A way to improve WRF-Noah and WRF-Alpine3D performances would be to increase the horizontal spatial resolution of WRF simulations (e.g., from 3 km down to 1 km), in order to force the snow cover models with more resolved atmospheric data.Increasing the WRF spatial resolution would also justify the activation in Alpine3D of the modules that take into account the aeolian transport of snow, that redistributes the snow between adjacent cells of the domain according to the forecasted wind field, as well as the impact of the local topography on the energy balance through the effects of incoming shortwave radiation shading and long wave radiation increase.The WRF-Alpine3D simulation could also be improved directly assimilating weather radar data in terms of snowfall rate over the study domain in order to provide more realistic snow accummulation patterns.The approach for weather data assimilation, in terms of point-wise nudging or statistical variational techniques, is an interesting open issue.https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.https://doi.org/10.5194/tc-2021-285Preprint.Discussion started: 3 November 2021 c Author(s) 2021.CC BY 4.0 License.

Figure 4 .
Figure 4. Flowchart of the model chain realized in this study.NCEP reanalyses are used to initialize the WRF model, which provides the atmospheric forcing data to drive WRF-Noah (online coupled) and WRF-Alpine3D (offline coupled) numerical models.

Figure 5 .
Figure 5.Comparison of in situ observed data and WRF simulations of daily air temperature (a), relative humidity (b), wind speed (c), incoming shortwave radiation (d) and precipitation (e).The shaded color of the dots indicate the altitude of the measuring sites, whereas the red line represents the best fit of the data and the dashed line shown indicate the bisector optimal line.

Figure 6 .
Figure 6.Time series of observed and simulated snow height (a) and snow height variation (b), averaged over all available stations.Observations are indicated by a solid gray line, while simulations are indicated by solid blue and red lines, corresponding to WRF-Noah and WRF-Alpine3D, respectively.The shaded areas correspond to an interval of ±1 standard deviation.

Figure 7 .
Figure 7.Comparison of simulated and measured daily means of snow height (a) and daily snow height variation (b), with blue and red dots indicating WRF-Noah and WRF-Alpine3D, respectively.The blue and red lines represent the best linear fit of WRF-Noah and WRF-Alpine3D data.

Figure 8 .
Figure 8. a) Observed and modelled frequency distribution of daily snow height variation at the 13 sites used in this study.The grey column represents the observed variations and the blue and red columns represent the snow height variations simulated with WRF-Noah and WRF-Alpine3D, respectively.b) ETS index for different snow thresholds with blue and red lines indicating WRF-Noah and WRF-Alpine3D scores, respectively.

Figure 9 .
Figure 9. a), b and c): MODIS, WRF-Noah and WRF-Alpine3D snow cover fraction maps.d), e) and f): MODIS, WRF-Noah and WRF-Alpine3D snow cover area maps derived by applying a threshold of 51% to the corresponding snow cover fraction maps.

Figure 10 .
Figure 10.CORINE broad leaved forest, coniferous forest, and mixed forest classes aggregated and reprojected from original 100 m resolution to 3 km resolution.

Figure 11 .
Figure 11.a), b and c): MODIS, WRF-Noah and WRF-Alpine3D snow cover duration maps derived from snow cover area maps.Fig. d) shows the differences of snow cover duration between WRF-Alpine3D and WRF-Noah, whereas Figs.e) and f) show the differences with MODIS of WRF-Noah and WRF-Alpine3D, respectively.

Figure 12 .
Figure 12. a) Time series of snow cover area fraction derived from MODIS (grey line), WRF-Noah (blue line) and WRF-Alpine3D (red line).b) Comparison of simulated and observed snow cover area fraction with blue and red dots indicating WRF-Noah and WRF-Alpine3D, respectively.The blue and red lines represent the best linear fit of WRF-Noah and WRF-Alpine3D data.

Figure 13 .
Figure 13.Time series of Jaccard index (a) and Average Symmetric Surface Distance index (b) for WRF-Noah (red line) and WRF-Alpine3D (blue line).

Figure 14 .
Figure 14.a) and b): WRF-Noah mean and maximum snow height.c) and d): WRF-Noah mean and maximum snow water equivalent.e) and f): WRF-Alpine3D mean and maximum snow height.g) and h): WRF-Alpine3D mean and maximum snow water equivalent.i) and j): WRF-Alpine3D and WRF-Noah differences of mean and maximum snow height.k) and l): WRF-Alpine3D and WRF-Noah differences of mean and maximum snow water equivalent.

Figure 15 .
Figure 15.Snow height (a) and snow water equivalent (b) averaged over simulation domain at low (< 800 m, blue lines), mid (800-1600 m, green lines) and high (≥ 1600 m, brown lines) elevation.WRF-Noah and WRF-Alpine3D predicted data are marked as dashed and continuous lines, respectively.

Table 2 .
Coordinates, real elevation and model elevation of the snow stations.

Table 3 .
Contingency table for the calculation and description of the ETS index.OY is observed yes, ON is observed no, FY is forecast yes,

Table 5 .
WRF-Noah and WRF-Alpine3D mean bias error (MBE), mean absolute error (MAE) and correlation coefficient (R) for simulated snow height and daily snow height variation.

Table 7 .
Mean values of Jaccard index and Average Symmetric Surface Distance for WRF-Noah and WRF-Alpine3D.