On the role of local and large-scale atmospheric variability in snow cover duration: a case study of Montevergine Observatory (Southern Italy)

Snow cover plays an important role in Earth’s climate, hydrological and biological systems as well as in socio-economical dynamics, especially in mountain regions. The objective of this work is to provide the first evidence about snow cover variability in the Italian Southern Apennines and investigate the forcing mechanisms controlling it. To this purpose, we present a new historical long-term (from 1931 to 2008) series of snow cover duration data observed at Montevergine Observatory, a mountainous site located at 1280 m above sea level. From the analysis of this series, it emerged a strong interannual variability, an overall reduction over time of snow cover days until mid-1990s and a recovery in the last 10-years. We model snow cover duration employing a multiple linear regression, considering both local and large-scale climate factors as explanatory variables. Our findings show that snow cover duration appears to be primarily dependent on temperature, which exhibits a positive trend in the considered time interval. However, the interannual and decadal fluctuations of the examined parameter are also strongly modulated by two large-scale patterns, the Arctic Oscillation and the Eastern Mediterranean Pattern. In the last segment of the considered time interval, the increase in temperature is not consistent with the dominant patterns of large-scale indices, which proved to be more effective in capturing the recent rebound in snow cover duration. The results demonstrate that snow cover duration is linked to the global warming by a non-trivial relationship and that its behaviour, in specific periods, can be largely independent from rising temperature tendency, according to the prevailing phase of large-scale atmospheric patterns.


Introduction
It is well known that snow cover is a key variable for several fields and applications and recently a great deal of attention has been paid to investigating the link between its variability and global warming. Snow cover is of crucial importance for climate studies, as it strongly affects the energy balance of a region by increasing land surface albedo (Sandells and Flocco 2014). It is a relevant input parameter for forecasting models in hydrology and water management (Terzago et al 2022). Moreover, its variations affect the mountain soil microbial processes as well as ecology and biogeochemical cycles (Freppaz et al 2012, Magnani et al 2017, Rogora et al 2018. Lastly, snow cover plays a fundamental role in agriculture, winter tourism and hydropower generation (Damm et al 2017, Wagner et al 2017. Thus, investigating the variations over time in the number of days with snow cover on the ground (i.e. snow cover duration), as well as of the main factors that explain its variability, is essential for determining potential future changes in climate that may heavily impact upon societies and economies (Beniston 2010). During the recent decades, satellite observations have opened up new opportunities for monitoring mountain snow cover (e.g. Bormann et al 2018, Tang et al 2022. However, they do not have yet a sufficient temporal extension for assessing long-term trends and variability. In this regard, despite their Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
well-known limitations (e.g. Notarnicola 2022), the in situ historical observations are still a cornerstone for climatological studies that search for evidence of past and recent changes in snow cover. However, as reported by Valt and Cianfarra (2010) and Diodato et al (2020), available historical long-term series on snow cover duration (hereafter, SCD) are limited worldwide (e.g. Hantel et al 2000, Dahe et al 2006, Brown 2019, Notarnicola 2022. Focusing on the Mediterranean region, which is one of the main climate change hotspots, most of the studies concerning on this variable are referred to the Alpine (e.g. Valt andCianfarra 2010, Matiu et al 2021) and Pyrenean regions (e.g. López-Moreno et al 2020, Bonsoms et al 2021). Instead, very few studies have been devoted to the Italian Apennines (Petriccione and Bricca 2019, Raparelli et al 2021, Diodato et al 2022. From these works, it emerged an overall negative trendency in SCD in almost all examined sites, in line with worldwide trends. However, the observed remarkable interannual and multi-year variability raised the need to improve the knowledge of the main factors (local and large-scale) driving snow cover dynamics (Notarnicola 2022).
In this study, we examine a new daily historical dataset of SCD measurements collected from 1931 to 2008 at Montevergine Observatory, in the Southern Italy. After discussing the evidences provided by this unique time series, we investigate the main factors controlling SCD variability through an empirical reduced-complexity statistical model. This paper is structured as follows: section 2 introduces the dataset and the statistical methods, section 3 shows the results obtained and in section 4 the conclusions are drawn.

Data and methods
2.1. Data 2.1.1. In-situ measurements The historical snow observations analysed in this study have been manually collected in the Montevergine Observatory (40.936502°N, 14.729150°E), a mountain site located at 1280 m above sea level in Southern Apennines in Italy. Specifically, snow depth (SD) measurements were collected on daily basis during the period from December 1931 to November 2008, with a few interruptions mainly located in the 1960s (Capozzi et al 2020). Unfortunately, the Observatory has no longer been manned since 2008 and the collection of several variables, including the snow depth, has been interrupted. Recent studies have employed snowfall data collected at Montevergine Observatory. For instance, Capozzi et al (2022) examined the winter variability of fresh snow depth for the period 1884-2020 and performed a cluster analysis to identify the meteorological scenarios that triggered the snow events. In Diodato et al (2022), annual snow cover duration values recorded at this site from 1952 to 1996 were used to develop a semi-empirical model relating this variable to meteorological conditions and elevation inputs.
In this paper, we examine SD in terms of duration, through the SCD index, obtained by summing the total number of days with snow cover on the ground (i.e. with SD > = 1.0 cm), for each year and for a given season. In the same period, homogenized 2m temperature (T) and total precipitation (P) observations collected in Montevergine Observatory are also available (Capozzi and Budillon 2017, Capozzi et al 2020).

Teleconnection patterns
To investigate about the relationship between large-scale atmospheric circulation and SCD variability in Montevergine, we have selected several large-scale patterns that are known to exert a relevant influence on Central Mediterranean area (Brunetti and Kutiel 2011, Caloiero et al 2011, Capozzi et al 2022. The large-scale patterns considered in this work are: the North Atlantic Oscillation (NAO), the Arctic Oscillation (AO), the Eastern Mediterranean Pattern (EMP), the Atlantic Multidecadal Variability (AMV), the Mediterranean Circulation (MC), the Western European Zonal Circulation (WEZC) and the North Sea-Caspian Pattern (NCP). The corresponding indices for these large-scale atmospheric circulation patterns have been retrieved from Climate Prediction Center of the NOAA's National Weather Service (Climate Prediction Center 2021). Each index is available on monthly basis and is standardized with respect to 1981-2010 climatology.

Methods
First, we extracted the main temporal variation modes of monthly SCD series using the singular spectrum analysis (SSA) method. This technique is a very useful tool that can be used for isolating periodic components and trends by decomposing the original series into the sum of a small number of independent and interpretable components (e.g. Rasouli et al 2020).
Subsequently, the available data (i.e. daily values of SCD, temperature, precipitation, and monthly largescale indices) have been aggregated according to different snow seasons: early season, defined as the period from November to January (NDJ), late season, from February to April (FMA), winter season, from December to February (DJF), full season, from November to April (NDJFMA). For temperature, precipitation and large-scale indices, seasonal values were obtained by simply averaging, for a given year, the corresponding monthly values; seasonal SCD values were computed from the daily values as a sum over different months.
A standard multiple linear regression analysis was carried out to quantify to what extent seasonal SCD variability can be explained by local variables (e.g. local T and P) and large-scale atmospheric patterns (Scherrer et al 2004).

Regression model
In order to simulate the behaviour in time of SCD observed in Montevergine, we have used an ordinary multiple linear regression model. In our case, the outcome variable is SCD anomaly (deviation from 1981-2010 mean), hereafter SCD a . To separately investigate the role of local and large-scale atmospheric variability, we have developed two multiple linear regression models: local and global. In the former, the explanatory variables are the temperature anomaly (T a ) and total precipitation amount anomaly (P a ) detected in Montevergine Observatory with respect to 1981-2010 baseline period. In the latter, we determined the most influential teleconnection indices using a strategy structured in the following two steps. Firstly, we calculated the Pearson correlation coefficient between SCD a and individual indices, and we selected indices more correlated with SCD a and less correlated with each other, to minimize multicollinearity issues. Then, we employ stepwise AIC selection to identify the final subset of predictors, obtained in the best performing specification (Yamashita et al 2007). We checked the adequacy of the resulting models by investigating potential violations of regression modelling assumptions (Poole andO'Farrell 1971, Dalic andTerzic 2021). The F-test was used to assess the overall statistical significance of the relationship (with p-value 0.05). To check the significance of individual regression coefficients, we adopted the t-test (Rubinfeld 2000). Finally, we investigated the individual contribution of the predictors to the model by various metrics of relative importance: Lindeman, Merenda and Gold (LMG) method, Beta, Last, First (Grömping 2007). We repeated the statistical procedure described above for all the considered seasons.

Model performance
In order to evaluate the modelling performance and measure the improvement in accuracy of global model compared with the local one, both the Mean Absolute Error (MAE) and the Mean Squared Error Skill Score (MSESS) were employed (Murphy 1988, Rust et al 2015: where SCD a k and SCD m k are the observed and model-predicted (local or global) values of snow cover duration anomaly at year k, respectively, n is the total number of available records, and MSE g and MSE l are the Mean Square Error (MSE) of global and local estimates, respectively. The MAE index measures the difference between the observed value and the outcome of each model, separately. MSESS is the proportion of improvement in accuracy of a prediction compared to a reference procedure, using MSE as the accuracy measure (Sutherland et al 2004).

Results and discussion
The first goal of this paper is to provide a characterization of both the short and long-term variability modes of SCD. Figure 1 shows the decomposition results obtained using the SSA technique. More specifically, the panels in figure 1 represent (starting from the upper) the original data, the long-term and decadal fluctuations, the seasonal variability mode, and the residual components of SCD time series for the 1932-2008 period (note that for convenience of calculation, we exclude the data collected in December 1931 from SSA and we remove the gaps in the time series). First, we see from figure 1(b) a decreasing long-term tendency of SCD signal and oscillations of a 10-year period. It is interesting to note that the amplitude of decadal variations gradually decrease until the late 1990s. As highlighted in panels (c) and (d) of figure 1, SCD is also subject to annual and sub-annual variations, related to regional climatic conditions, and reveals a very high interannual and intermonthly variability typical of the mid-latitudes. A detailed quantitative analysis of Montevergine SCD time series confirms a general reduction in SCD during the investigated period for each snow season. In particular, SCD decreased by 12%, 17%, 12%, and 26% in the period 1969/70-2007/08 compared to the previous 30-year period (1932/33-1963/64), in full, early, late and winter season, respectively. The substantial decline in SCD is consistent with previous studies (e.g. Valt andCianfarra 2010, Klein et al 2016). It is worth noting that the most pronounced reduction in SCD occurs in the winter season (figure 2(a)), when the most relevant increase in T is observed ( figure 2(b)). In fact, although the T time series measured at Montevergine Observatory exhibits a positive trend in all seasons, it is statistically significant (at 95% confidence level) only in late and winter season (+0.12°C decade −1 and +0.18°C decade −1 , respectively). This result is in line with some previous studies, even   [1998][1999][2000][2001][2002][2003][2004][2005][2006][2007], despite a further slight increase in the average T. More specifically, since 1998 the SCD has rebounded to values comparable to those observed in previous periods (1948-1957, 1958-1967, 1978-1987) where the average T was about 1.0°C lower than 1998-2007. It is worth noting that a recovery (generally modest) of SCD or some other important snow climate indicators from the late 1990s or early 2000s has been also observed in several other locations on European Alps (Scherrer et al 2013, Marcolini et al 2017, Matiu et al 2021. Following such evidence, we conducted a multiple linear regression analysis to better understand the potential main drivers of snow cover variability in Montevergine. Table 1 presents the Pearson correlation coefficients between all considered predictors and SCD a , whereas table 2 shows the coefficients of final multiple regression models and their significance level (according to t-test results described in section 2). Given that all assumptions are sufficiently fulfilled for both local and global models in all seasons, we conclude that a meaningful linear relationship between SCD a and the predictors exists.
In particular, the proportion of total variation in SCD a explained by T a and P a in local model ranges from 40% (in full and early seasons) to 57% (in winter season), as shown in table 3. It is worth observing that T a provides a larger relative contribution compared to P a in all models, although the percentage varies slightly across seasons. In fact, for all metrics and snow seasons, T a relative importance is around 80% of the total R 2 . This result is in line with previous studies, in which the temperature has been identified as the main driver (with respect to precipitation) influencing snow cover in regions where the mean winter T is not far below the zerodegree limit (Dahe et al 2006, Brown 2019, Diodato et al 2022. Regarding the global model, Pearson's pairwise correlation analysis between SCD a and large-scale indices (table 1) shows the presence of statistically significant positive and negative correlations with EMP and AO indices, respectively, in all seasons (except in early season for AO index). These indices appeared to be the main significant predictors for SCD a ; in winter, NCP index also turned out significant (table 2). Overall, the explanatory variables in the global model account for no more than 50% SCD a variance (table 3). More in detail, EMP index is the most influent (in terms of relative importance metrics introduced in section 2) in full, early and winter season, while in the late season, a majority (around 60%) of the total explained variance is accounted by the AO index. It is useful bearing in mind that the negative AO phase (associated with a weak polar vortex allowing cold air outbreaks across mid-latitudes (Thompson and Wallace 1998)) and the positive EMP pattern (characterized by positive 500-hPa geopotential height anomalies over northern Atlantic and by negative anomalies over Central and Eastern Mediterranean basins (Hatzaki et al (2007)) are both synonymous of synoptic weather types that promote snowfall events in Montevergine as recently demonstrated in Capozzi et al (2022). The strong impact that AO and EMP indices have on local weather conditions is highlighted in table 4. It is evident that the EMP index is well correlated with temperature in all seasons, especially in winter, whereas the AO index has a good linkage with the precipitation. In summary, negative AO and positive EMP regimes determine the 'ideal' local weather conditions for snow occurrence and persistence in the southern Apennines, whereas in the opposite AO and EMP phases (positive and negative, respectively) the 'worst' conditions are generally observed. However, it is necessary to consider that temperature and precipitation variability observed in Montevergine is conditioned not only by the large-scale atmospheric patterns (such as EMP and AO), which act on different time scales (sub-monthly, seasonal, interannual and decadal), but also by strictly local factors (i.e. site exposure, orography, etc), and, on longer periods, by global warming. Therefore, the two linear regression models (local and global) provide different information for our analysis concerning the forcing of SCD variability.
As revealed by figures 2(c) and (d), from early 1930s to the end of 1960s, negative AO phases and positive EMP prevailed in winter season. In the subsequent 20-years period, a positive (negative) tendency of AO (EMP) occurred, culminating between late 1980s and early 1990s, when both indices were in a phase unfavourable to snow events in Montevergine. Finally, a reverse trend for both indices has been observed in the last time interval.
The comparison between local and global models shows that the former generally performs slightly better than the latter in terms of both MAE and MSESS (by 8% and 0.2 on seasonal average, respectively).
More compelling evidence about the performance of the models is provided by figure 3, which shows the comparison between the estimated and the observed SCD a series during the period from 1931 to 2008 in full season ( figure 3(a)), early season ( figure 3(b)), winter season (figure 3(c)) and late season ( figure 3(d)). In order to remove the high-frequency variability, both simulated and measured signals have been smoothed with a 10years lowess filter. It is worth noting that both local (marked as red solid line) and global (highlighted as dashed magenta line) models are able to replicate the time pattern of the observed SCD a (blue solid line) in an acceptable way, especially in full ( figure 3(a)) and winter (figure 3(c)) seasons.
However, a detailed inspection of figure 3 reveals a considerable difference in model performances in the last segment of the analysed time interval. More specifically, the global model reproduces better than local one the rebound of SCD observed from late 1990s. To quantify this difference, we have computed the MAE and MSESS for separate sub-periods, having a length of 10 years (except for the first one, which spans over 7 years). The results, sketched in table 5, show that in the 1998-2007 period the global model outperforms the local one in full, early and winter seasons, as shown by the lower MAE values and by the positive MSESS score. On the other hand, the two models display a similar performance in the late season. It is also interesting to note that in the first three sub-periods the local model generally outperforms the global one in all the considered seasons. Subsequently, in full and winter seasons the model performances exhibit a decadal behaviour, well synthetized by the regular switching of MSESS score from positive to negative values. The early season is aligned to this temporal pattern In late season, the local model performs generally better than global one, except in the 1988-97 period, in which the large-scale teleconnection patterns are more effective in reproducing the SCD a variability (in contrast with the other seasons). In full, winter and early seasons, the behaviour of the models appears to be strictly dependent on the relationship between the considered local and large-scale atmospheric variables. From early 1930s until the late 1960s, the synoptic atmospheric circulation was very favourable for snowpack persistence (negative AO and positive EMP). In this period, the observed negative T and positive SCD anomalies are consistent with large-scale patterns. In the last 40-years of the considered time interval, the gradual increase in T and decrease in SCD is again consistent with the positive (negative) trend of AO (EMP) until early 1990s. Subsequently, the rising temperature tendency was not consistent with the dominant neutral AO and positive EMP patterns, which determined a recovery in SCD. In other words, in the last 15 years of the investigated period, the increase in temperature induced by the global warming can be clearly distinguished from the large-scale atmospheric  circulation variability. Therefore, the local regression model fails in correctly reproducing the recovery in snow cover duration because it is trained on seasonal temperature anomalies, which, due to the global warming, are mostly above the climatological average. However, it is important to note that the worse performance of the local model may also be ascribed to sub-seasonal variations in temperature and precipitation (induced by the largescale atmospheric circulation) that are not taken into account by the model. In the analysed period , the global warming signature in Montevergine data is very strong, especially in winter season, as previously discussed. In the same period, for AO and EMP indices we found a negligible positive (+0.10 decade −1 ) and negative (−0.02 decade −1 ) tendency, respectively. Thus, in the considered period there is a statistical independence between large-scale indices (AO and EMP) and temperature trend, which does not imply complete physical independence (Kim et al 2022). On the basis of this assumption, it is possible to separate the contribution of large-scale indices and global warming on snow cover duration variability.
From a physical point of view, the AO and EMP patterns may be considered as 'internal' modes of climate system, whose long-term and/or decadal variability may be not necessarily aligned to the global warming trend (and, therefore to local temperature anomalies).
This aspect has been previously highlighted in a broader context by Cohen and Barlow (2005), who concluded that the positive temperature trend recently observed in northern hemisphere is largely independent of the AO and NAO indices.
According to the findings of our study, on decadal scale the behaviour in time of SCD might take directions that strongly diverge from temperature evolution and that are closely connected with specific large-scale atmospheric patterns. The latter are notoriously characterized by a very complex variability over the time, which may amplify the global warming effects (as happened, according to this study, in the 1970-1990 period) or strongly mitigate them, as we observed from late 1990s to late 2000s.

Conclusion
Snow cover is a crucial variable in several fields and has been officially defined an Essential Climate Variable by the Global Climate Observing System (WMO 2010). This work aims to provide new insights about past snow cover variability, analysing a 77-year long series of snow cover duration measurements recorded at Montevergine Observatory in the Apennines, an area that has been very scarcely investigated so far. Our results show an overall decrease over time of snow cover duration, most notable in winter season, until the mid of 1990s and then a rebound in the last 10-years period, as well as strong decadal and interannual variabilities. We developed two empirical multiple linear regression models that allow disentangling the main factors affecting SCD variability. Our results show that long-term reduction in snow cover duration at Montevergine can be mainly ascribed to the gradual increase in temperature. However, on shorter time scales we show that the SCD fluctuations are also strongly modulated by AO and EMP teleconnection patterns. In the last 10 years of our sample data, the synoptic scale circulation patterns combined for conditions favourable to extended period of snow cover despite the background warming trend.
The findings of our study indicate, in line with Brown (2019), that occasional periods of longer SCD, related to main annual or decadal mode of teleconnection indices, might be expected despite projections of long-term declines of SCD due to global warming.
Starting from those results, future research activities should be primarily addressed to investigate about the physical mechanism connecting the rising temperature tendency and the large-scale atmospheric variability driving the SCD variability. This aspect is essential to improve the performance of climate models designed to delineate future SCD evolution. Further efforts should be also devoted, on one hand, to replicate our approach in other mountainous regions of Mediterranean area (if possible including data for the last 15 years). On the other hand, in order to increase confidence on our results, in future studies we will evaluate empirical models of increased complexity to reconstruct the SCD variability.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.