# Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data

Nazzareno Diodato, Sara Rubinetti, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Sara Rubinetti, Gianni Bellocchi. Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data. Hydrological Sciences Journal, 2024, 70 (2), pp.294-310. 10.1080/02626667.2024.2427356. hal-04954567

HAL Id: hal-04954567

https://hal.inrae.fr/hal-04954567v1

Submitted on 18 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data

Nazzareno Diodato, Sara Rubinetti & Gianni Bellocchi

**To cite this article:** Nazzareno Diodato, Sara Rubinetti & Gianni Bellocchi (07 Nov 2024): Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data, Hydrological Sciences Journal, DOI: 10.1080/02626667.2024.2427356

**To link to this article:** https://doi.org/10.1080/02626667.2024.2427356

<figure><img src="image_2.png" /></figure>

Accepted author version posted online: 07 Nov 2024.

<figure><img src="image_3.png" /></figure>

Submit your article to this journal

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

View related articles

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

CrossMark

View Crossmark data

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /></figure>

**Publisher:** Taylor & Francis & IAHS

**Journal:** *Hydrological Sciences Journal*

**DOI:** 10.1080/02626667.2024.2427356

# Enhancing decadal snowfall forecasts in the Mediterranean mountains through informed atmospheric variability and climate data

Nazzareno Diodato<sup>a</sup>, Sara Rubinetti<sup>b</sup>, Gianni Bellocchi<sup>*a,c</sup>

<sup>a</sup> Met European Research Observatory – International Affiliates Program of the University Corporation for Atmospheric Research, Via del Pino 47, 82100 Benevento, Italy. ORCID: https://orcid.org/0000-0001-9549-1583

<sup>b</sup> Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, List/Sylt, Germany. ORCID: https://orcid.org/0000-0002-2236-6113

<sup>a,c</sup> Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France. ORCID: https://orcid.org/0000-0003-2712-7979

*Corresponding Author:* Gianni Bellocchi, Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France. E-mail: gianni.bellocchi@inrae.fr
<figure><img src="image_10.png" /></figure>
**Abstract**

Forecasting decadal-scale snowfall, crucial for global water management, is challenging due to the complex interplay of environmental factors. This study projects the number of snowfall days (NSD) up to 2060, using the extensive time-series data from the Montevergine Observatory, southern Italy (1884-2022). We pioneered an innovative statistical model, accounting for antecedent time-lags, and exogenous support from the Teleconnection-Climate Pattern Index, incorporating large-scale (Arctic Oscillation) and smaller-scale (temperature) forcings. Our projections reveal the influence of decadal and multidecadal oscillations throughout the forecast period and suggest an increase in NSD after 2030, notably shifting in the 2040s to 2050s, averaging from about 20 to 30 snowfall days annually. The frequency of snowfall deficit years (reaching -1 standard deviation) remains, however, high in the first part of the forecast. Despite the limitation of a single-site study, this trend is consistent with projections from various regional circulation models for increased extreme snowloads in Italy.

**Keywords:** Arctic oscillation; autoregressive model; decadal predictions; snowfall

**1. Introduction**

Decadal climate prediction is widely recognised as a key planning horizon for societal entities, benefiting from various types of information (Solaraju-Murali *et al.* 2021). However, long-term forecasts remain a significant challenge due to the predictability of variables such as sea-surface and land temperatures (Sandgathe *et al.* 2020), and precipitation, in particular concerning the snowfall component (Henderson *et al.* 2018). Precipitation prediction is complicated also by the orography (Schroeer and Kirchengast 2018; Smith and Bookhagen
2018) and intricate modelling of the interactions among different spatial scales of atmospheric
modes affecting precipitation (Efthymiadis et al. 2007) and climate evolution (Beniston et al.
2018). The complex regional variability in climate affects both snowfall and snow cover by
altering how water resources are partitioned in the flow pathways of the global water cycle
(Hale et al. 2023; Figure 1). Thus, snowfall and snow depth play a crucial role in shaping the
interannual variability of water resources and the hydrological cycle (Wang et al. 2018),
particularly in regions where water storage capacities are insufficient (Barnett et al. 2005) and
in elevated catchments, which are more sensitive to environmental water storage (Beniston et
al., 2014; Marty et al. 2017; Musselman et al. 2021). It is essential to dedicate more effort to
understanding the climatic properties of various mountain regions and the links between
global climate forcing and regional climate responses (Beniston 1997).

Historical time-series assessments indicate a decline in snow cover (Klein et al. 2016), snow depth (Meredith et al. 2019) and number of snowfall days in the European pre-Alpine region (Diodato et al. 2020), as well as at different latitudes worldwide (Shi and Liu 2021). Concerning the central Mediterranean area, the change in snow response to climatic oscillations can be traced back to end of the Little Ice Age (LIA) (Diodato et al. 2019; Figure 2a).

Many studies on snowfall variability (Dye 2002) have demonstrated correlations with local meteorological variables (e.g., precipitation, temperature, geopotential height) and large-scale circulation influences (Bednorz 2004). For example, Notarnicola (2020) found a negative correlation between temperature and snow cover in the winter season in Eastern Europe. Kis and Pongrácz (2021) discovered a significant correlation between temperature and snowfall (p<0.05) in all European areas from October to May. Recent research has expanded the focus on the number of snowfall days (NSD, defined as the day with snowfall ≥ 0.01 m) and climate change to include the impact of teleconnections and large-scale climatic events on snowfall.
Bailey *et al.* (2021) examined European large-scale circulation anomalies, associating phases of negative Arctic Oscillation (AO) with cold air advection and heavy snowfall in mid-latitudes. The AO, a mode of climatic variability, influences precipitation and temperature in mid-latitudes (Hamouda *et al.* 2021).

Capozzi *et al.* (2022) detected a notable decrease in snowfall frequency between the 1970s and 1990s for Montevergine, mainly attributed to a strong positive trend in AO. Additionally, recent analysis by Yang *et al.* (2023) observed that the negative phase of AO leads to a deepened East Asian trough and a small trough over Europe, resulting in colder-than-normal Eurasian mid-latitude weather and favourable conditions for snowfall events.

The number of snowfall days contributes to a significant source of change, prompting inquiries into whether trends in atmospheric and climatic factors drive snowfall trends. It raises the question of whether historical snowfall data can effectively inform future projections. The Greek historian Polybius (c. 206 – c. 118 BC), in his Histories, emphasises the enduring relevance of understanding past patterns for anticipating future developments (Martines 2018):

> Thanks to the memory of the past, [...] evaluate any situation with greater certainty because the knowledge of what has already happened makes predictions for the future more solid and well-founded.

Deterministic forecasts derived from empirical (time-series based) or dynamical models can offer potential solutions, contingent upon how effectively the model aligns with its intended purposes, and serves as a practical guide for various applications. Physically-based models, which are based on the numerical representation of the key physical and chemical processes, while capable of tracing water movements to anticipate spatiotemporal trends, are cautioned against due to their inherent complexity and high processing costs. Actually,
the assessment, in probabilistic sense, of regional climate changes requires large ensemble
simulations (Xie et al. 2015). Global Circulation Models (GCMs) and Regional Circulation
Models (RCMs), despite their development and enhancement for broader climate
projections, introduce uncertainty, as illustrated by Nerantzaki et al. (2020) in their work on
the combined uncertainty of five RCMs and a hydrological prediction framework. Yuan et
al. (2023) emphasised the significant role of GCMs and RCMs in hydrological projections but
acknowledge the inherent uncertainty associated to the internal variability arising from the
intrinsic chaotic nature of climate. Consequently, statistical forecasts based on local
information could serve as a valuable benchmark to assess the added value of dynamical
forecasts, providing a complementary perspective (Tommasi et al. 2017). Simpler statistical
approaches, in particular, may offer supplementary insights into how the future state of a
system is influenced by past NSD information and potential predictors. These approaches
implicitly examine all factors impacting hydrological data evolution, revealing complex
non-linear interactions between input and output variables (Vishwakarma et al. 2023).
Although the hydrological study of snow processes poses challenges due to the complex and
varied orography and terrain aspects, the temporal trend of snowfall across relief strongly
responds to climatic fluctuations, especially in precipitation (Misra et al. 2020).

Research on precipitation patterns in mountainous locations is pertinent due to their
untamed and frequently uncontaminated natural ecosystems (Capozzi et al. 2022). Merino et
al. (2014) focussed on the Iberian Peninsula, specifically the Castilla y Leon plateau,
highlighting a negative trend in NSD from 1960 to 2011 and revealing interesting links
between regional snowfall variability and large-scale atmospheric patterns. Baltaci et al.
(2020) observed positive trends in NSD in the central Anatolian plateau (Turkey) from 1967
to 2006. Conversely, Blahušiaková et al. (2020) recorded negative changes in snowfall days
in mountainous catchments in Central Europe between 1966 and 2012.
Brázdil et al. (2023) confirmed similar patterns for snowfall days during winters (December–February) averaged for the territory of the Czech Republic from 1961 to 2021. Most of the snow extents in southern Europe have been consistently decreasing since the beginning of the third millennium (Figure 2b). Their extension remain confined to continental regions, including east France, the Balkan peninsula and Turkey (Figure 2c). Some small areas of snow cover persist throughout southern Europe, particularly in the Alps at higher elevations and, secondarily, the Apennines Mountains.

Focusing on the Italian peninsula, the distribution and hydrological importance of snow-cover duration (Valt and Cianfarra 2010; Enzi et al. 2014; Fratianni et al. 2015; Diodato et al. 2019) and snow depth (Nigrelli and Collimedaglia 2012; Terzago et al. 2013; Marcolini et al. 2017) are extensively studied over multidecadal and longer periods, whereas NSD has received comparatively less attention.

The study of time-series (Time Series Analysis, TSA) has demonstrated considerable success in reproducing observed climate variability. Different works employing TSA approaches in forecasting hydrological time-series have shown the utility of Autoregressive Integrated Moving Average (ARIMA) models, which remain widely used in TSA (Box et al. 2013). Refined approaches to time-series forecasting techniques have emerged in recent years. Montgomery et al. (2015) suggested that by incorporating periodic components, seasonal ARIMA (SARIMA) models can be created. Rubinetti et al. (2020) provided a two-decadal forecast of the Po River's discharge in Northern Italy, utilising a statistical technique that combines the independent implementation and cross-validation of autoregressive models and neural networks preceded by the denoising of the time-series.

This study introduces a pioneering approach by employing TSA with novel exogenous (X) inputs to enhance decadal predictions in a hydrological context. In this way, we have developed a time-varying (TVAR) Periodic AutoRegressive Moving Average (PARMA) process (referred
to in line with Wang *et al.* 2005) with Exponential Generalised Autoregressive Conditional Heteroscedasticity (EGARCH) - PARMAX(TVAR)-EGARCH - to predict NSD for the longest observed Mediterranean time-series (1884-2022) at the Montevergine Observatory in the southern Apennines of Italy. While ARMA processes are commonly applied in hydrological studies, GARCH models are widely used in econometrics (e.g., Gao *et al.* 2018; Fuente *et al.* 2020). Nevertheless, their applicability in climate studies has already been explored (e.g., Wang *et al.* 2012; Diodato *et al.* 2023).

The performance of ARMA-based models can significantly improve when external forcings are introduced in the models through the use of exogeneous variables (e.g., Sarhadi *et al.* 2014). In this study, the model incorporated exogenous climatic inputs to account for both large-scale (Arctic Oscillation teleconnection) and smaller-scale (temperature) climate forcing agents, aiding to mitigate uncertainty in interannual and longer-term trend variability. The PARMAX(TVAR)-EGARCH model underwent development using data from the training period 1884-1991, validation for the interval 1992-2022 and application for forecasts in the time-horizon 2023-2060. The validation period was constrained to 31 years due to the inability to deduct additional years from the training period. This limitation was necessary to ensure that an adequate number of years remained to process the selected model type. The deliberate decision to extend the duration of the forecast period beyond that of the validation period was made to provide a broader time window for comparing the model's predictions with climate model estimates, which are often projected somewhere between the mid-21st century (around 2050) and the end of the century (around 2100) (e.g., Faggian *et al.* 2024).

## 2. Environmental setting and model building

### 2.1 Study area and snow climatology
The Montevergine mountain, also known as Partenio, is situated in the Campanian sub- Apennine range in southern Italy (Figure 3a, b). It occupies an eccentric position near the ridge chain that divides the Italian peninsula into the Adriatic and Tyrrhenian segments (Figure 3b, c). The Montevergine Observatory is located found 40° 56' North and 14° 43' East, at an elevation of 1270 m a.s.l., providing a view of the Sabato Valley (Figure 3d).

From the Observatory's vantage point, the panoramic view extends over 100 km, encompassing the regions of Molise and Apulia to the east, and Basilicata to the south. On exceptionally clear days, the view extends to the coasts of the Balkans. Nestled in the free troposphere and remote form human influence, the Montevergine Observatory experiences a Mediterranean climate characterised by strong seasonal variability. Rainfall is minimal in summer and reaches its peak in winter (November-March, Diodato et al. 2021), when precipitation takes on a snowy character. Figure 3b illustrates the annual distribution of snow depth in Italy averaged over the period 1921-1960, showing an increase toward the northwest, the Alps, and then over the peninsula along the Apenninic chain, where variations in snow depth are most pronounced.

In the recent decades a significant decline in snowfall occurred and, consequently, the historical snowfall values shown in Figure 3 may be overestimated; anyway, they still serve as a valuable documentary source for studying this phenomenon. Unlike the early 20th century, when numerous survey stations were established by the former SIMN to monitor daily accumulations, todays such surveys are conducted by only a handful of amateur observers. Consequently, constructing an updated map comparable to the one presented in Figure 3b is no longer feasible. The effect of exposure is evident on the slopes of the Apennine region. During the winter season, prevailing westerly and south-westerly winds on the Tyrrhenian side, and northern winds on the Adriatic side, result in a significant difference in snowfall. Snow is scarce on the Tyrrhenian side but more abundant on the Adriatic side. The moderating influence of the sea
significantly affects the coasts of central Italy, southern Apulia, areas around Sardinia, and a substantial portion of Sicily. In peninsular Italy, the most conspicuous snowfalls, though notably lower than those reported for the Alps, occur in an extensive area of the Central Apennines, encompassing the highest mountain formations across the entire Apennine chain. In this area, in fact, the isonivals of 400 cm encircle the highest peaks, generally exceeding altitudes of 2500 m a.s.l., including the Sibillini Mountains, Gran Sasso of Italy, Terminillo and other mountain groups. Particularly in the southern Apennine areas, snow values decrease further, with 200 cm are slightly exceeded only in Campania, on Monte Partenio (including Montevergine), and in Calabria, on Monte Pollino (Gazzolo 1973).

Figure 3c focusses on the distribution of the annual mean snowfall in the vicinity of the Montevergine site, featuring circular isonivals characteristic of Montevergine with a snow depth interval of 100-300 cm. The annual means of days with snowfall roughly align with the isonivals. In the central and southern Apennines, the frequency at the same altitude is notable, as winter in these areas, unlike in the Alps, constitutes the season with the highest rainfall of the year. Consequently, there are 30 days of snow at Terminillo (2217 m), 60 at the Campo Imperatore station (2125 m) on the Gran Sasso, and 39 at Montevergine (1270 m), in comparison to the Lago Cereser station, situated at 2500 m a.s.l. in the Western Alps, which records 47 days of snowfall (Pinna, 1973).

While snow may not typically play a prominent role in the water balance of these areas, it consistently occurs almost every year, with a mean annual NSD estimated between 35 and 40 d yr⁻¹ (Pinna 1973). Despite the relatively modest coverage of mean annual NSD, its significance arises in specific situations, such as impacting the spring discharge of river basins (Diodato et al., 2023). This importance is attributed to the fact that the duration of snow cover is influenced not only by average intensity but also on the frequency of snowfall (Pinna 1973).
## 2.2 Snowfall data collection

The MVOBS snowfall data encompasses daily observations spanning from 1884 to 2022. These records, maintained by successive Benedictine Fathers responsible for this service, reveal a legacy of meteorological dedication. Father Giuseppe Llobet, the inaugural observer from 1884 to 1919, directed the initial meteorological efforts at the Montevergine Observatory (Figure 3e). Notably, there is a gap in evidence before 1893 regarding the use of the NSD measurement. Beginning in 1894, snow data found its repository in the *Giardinetto dell'Ave Maria*, a cloister located close to the northern side of the meteorological tower (Capozzi and Budillon 2017; Capozzi *et al.* 2022). The custodial baton passed through the hands of Fathers Ildebrando Mancini, Ilario Mauro, Ugo Inzan and Giulio Corvino, who sustained observations from 1920 to 1938. In 1939, Father Virginio Cinella assumed responsibility for the Observatory's services, marking the onset of an intense period that ensured until 1961. Subsequently, a scaled-back service managed observations in the ensuing biennium, culminating in their discontinuation in 1964. Father Amato Gubitosa, the Observatory's director from 1968 to 2006, significantly contributed of this study by extracting data from the original registers. These records were processed using metadata documented in the Observatory's chronicles. From 2007 onwards, Father Benedetto Komar, the Observatory's director, oversaw renovations and the installation of new hardware, including digital instruments. Post-2006 updates draw from Capozzi *et al.* (2022) and the Montevergine Observatory website (http://www.mvobsv.org), as provided annually.

## 2.3 Explanatory data analysis and predictability

The exploratory data analysis began by examining the NSD time-series spanning from 1884 to 2022. This included scrutinising statistical data to identify irregular or random
fluctuations, as well as variations in mean and variance. Furthermore, an assessment of
fluctuation patterns, mean and variance changes, and the presence of heteroscedasticity is
crucial to demonstrate stationarity.

Figure 4a shows the progression of the annual NSD, with the Mann-Kendall test (indicated by the gray dashed linear trend) robustly supporting a significant decreasing trend (p<0.01). Figure 4b illustrates the NSD residuals, which exhibit volatility clustering and suggest a time-varying and significantly decreasing conditional variance (Figure 4c). However, the time-varying component included in the autoregressive model we applied overcomes this issue, as it is also suitable for time-series affected by trends.

The data exploration involves also visualising the autocorrelation function (ACF) and partial autocorrelation function (PACF), conducting a dynamic analysis in the space-time domain. In Figure 5a and 5b, the computed ACF illustrates the correlation between NSD and different delays over residual values. The lag k represents the correlation between the standardised residuals of NSD at time t and time t+k, aiding in determining the autoregressive model order needed for optimal data fitting. Additionally, 95% confidence limits assess the statistical significance: If the estimated coefficient falls outside the probability boundaries for a given lag, a statistically significant correlation at the 95% confidence level is indicated. Similarly, PACF illustrates the correlation at a certain k, after removing the effects of the previous lags. The ACF (Figure 5a) and PACF (Figure 5b) delineate the pattern of linearly decreasing autocorrelation in residuals and an exponential pattern, respectively. In the ACF, the first eight spikes surpass the critical boundary, while in the PACF, only the first two spikes exhibit a similar pattern.

According to Valenzuela *et al.* (2008), discernible patterns in the ACF and PACF plots indicate the potential suitability of an ARMA(P,Q) model with initial order parameters (P and Q) set to 1 from the initialisation component model. The model was first trained and then Q and P were arranged based on the validation results.
The examination of the residuals sequence (Figure 4b) involves the McLeod and Li (1983) test (ML). The test results indicate that ML has a *p*-value of 0.0148, which is below the specified control *p* (0.05). Moreover, the NSD time-series seems non-stationary in both first and second order. This prompted our investigation towards models where variance and autocorrelation do not remain constant over time, and heteroscedasticity could be incorporated in the ARCH model at the same time. This suggested that the time-series was heteroskedastic, which led us to explore the appropriate PARMAX(Q,P) × EGARCH(Q,P) process. The validation successively confirmed the model to fit the NSD equal to PARMAX(1,1)*m* × EGARCH(1,1), where *m*=28 represents the periodic cycle associated with the *location* (L) component. The detailed explanation of each component of the model are reported in Section 2.5. Cycle *m* associated with the *location* component of the model was initially searched for by spectral density function and refined during the validation stage. Spectral density serves as a metric of how a signal's power is distributed across various frequencies. While the spectral density and the ACF contain the same information, they present it in distinct manners: Spectral analysis facilitates the identification of inherent periodicities, enabling the detection of cycle components in a time-series, along with the lag associated with these periodic pattern (e.g., Azari *et al.* 2021).

Additional exploratory data analysis was conducted based on the Hurst (H) exponent, a crucial
metric for establishing the memory strength of a time-series and its forecasting potential. The
exponent, where 0≤H≤1, defines the rate of stochastic processes within a time-series. It is
intricately linked to the fractal dimension (D) of the time-series through the relation D=2−H,
providing insight into its autocorrelation or long-term memory properties. An H value greater than
0.5 signifies persistent or trending behavior in the time-series, akin to a random walk, indicating
long memory effects. Conversely, an H value of 0.5 indicates an anti-persistent time-series signal.
In such cases, if the time-series increased in the previous period, it is more likely to decrease in
the following period, rendering it more irregular than a purely random signal. To estimate H, we
used the SELFIS self-similarity analysis software (http://www.cs.ucr.edu/tkarag/Selfis/Selfis.html), and applied the non-parametric Rescaled Range Analysis method (Mandelbrot and Taqqu 1979) to the detrended NSD time-series. This method is known for its robust estimation capabilities in the presence of stable distribution (Taheriyoun and Moghimbeygi 2017).

To achieve this, we initially employed the source time-series to determine the Hurst exponent (H), aiming to discern the presence of a long memory effect, its nature, and whether this memory could extend into the future during forecasting. This computed exponent H for the Montevergine instance is 0.90. This indicates a positive outcome, signifying that if the time-series has been growing in the past, it is likely to continue doing so in the future, and vice versa. Moreover, as H exceeds 0.65, the exponent suggests that this memory can indeed be used for future expansion.

For a more in-depth analysis of the predicability, the attractor offers insights into the reliability of the forecast. An attractor is the set towards which a dynamically disturbed system evolves after an extended period. The reliability of the forecast hinges on whether the system evolves in an orderly manner following the disruption. Dependability criteria for the forecast of the NSD time-series were subsequently calculated, aivotal stage in TSA. In the time-state domain, the attractor for Montevergine snowfall days exhibit a more tightly wapped directional branching (Figure 5c, red ellipse), indicating a potentially predictable pattern, provided reasonably accurate projections are attainable.

## 2.4 PARMAX(TVAR)-EGARCH model architecture

The greatest advantages of using such autoregressive models become apparent when dealing with the complexities of predicting dynamics influenced by flow periodicity and lag time in non-stationary and unpredictable time-series, such as NSD. These models excel in situations
where the highest amounts of information transfer between states at time t-1 and time t is crucial. Figure 6 presents a visual representation of the forecasting methodology, outlining the key stages of the PARMAX(TVAR)-EGARCH modelling applied in this case study. The random variable y at time t (i.e., y_t) is modelled to depend not only on y_{t-i} (where i is 1) but also on a fixed period m in the initial phase. A specific technique for incorporating the periodic (m) component into the autoregressive model involves an additive structure, which is essentially composed of two components. The first is interpreted as the location component, as described by Montgomery et al. (2015):

$$y_t = S_{t-i}(\text{ARMAX}(1,1)) + m_{t-jm} \quad (1)$$

where y_t is the time-series of NSD to be forecasted at time t over a period m; S_{t-1} represents the random component, described as an ARMAX (autoregressive exogenous) process with orders p=1 and q=1, and m_t denotes the deterministic periodic cycle with period m. The parameter specifications for the location component are given by:

$$\mu(L) = \omega + \text{AR}_1(1) \cdot \Phi^{(t-1)} + \text{AR}_2(1) \cdot \Phi^{(t-2)} + \text{Periodic}(m) + X(t) \cdot \beta_L + \epsilon(t) \quad (2)$$

Here, ω (constant term) equal to 8.8519 (initial value 1.8368) represents the baseline value of the time-series when all other factors are zero. AR_1(1) and AR_2(1) are the first and second autoregressive components. Φ (autoregressive parameters) equal to 0.1722 (initial value 0.0200) and to 0.7759 (initial value 0.9500) control the decay of the influence of past values with increasing lag (i). Periodic(m) represents a seasonal component with a period m=28, X(t) is the exogenous variable with the associated parameter β_L, and ε(t) reflects the difference between the actual value at the time t and the models' predicted value, capturing the unexplained component after accounting for the effects of past values, seasonality and exogenous variables.

In score-driven models, the residual term ε(t) is often analysed through the score term, which acts similarly to the gradient of the log-likelihood function with respect to a time-
varying parameter. This score term guides dynamic adjustments to the parameters based on observed data, enabling the model to adapt and potentially deviate from the assumption of normality for the error term. The final score component is the score vector of the (predictive) density of the observed time-series, highlighting the model's ability to handle deviations from normal distribution for the irregular component (Lit et al. 2022). The score-driven model has three main advantages: (i) the 'filtered' estimates of the time-varying parameter are optimal in terms of Kullback-Leibler divergence, providing an optimal filtering mechanism with respect to the true data generating process; (ii) the likelihood of these observation-driven models is established in closed form; and (iii) the predictive ability is comparable to that of parameter-driven competitors (Koopman et al. 2016).

For the log-scale component, which accounts for heteroscedasticity in score-driven frameworks, the specification is:

$$ \sigma(S) = \exp(\omega + AR_2(1) \cdot \Phi^{(t-2)} + X(t) \cdot \beta_S + \epsilon(t)) \quad (3) $$

where $\sigma$ represents the standard deviation, $AR_2(1)\Phi^{(t-2)}$ is the second autoregressive component of the first order and $X(t)$ is the exogenous variable with the associated parameter $\beta_S$. The residual term $\epsilon(t)$ in this equation has the same score-driven meaning as in the location component, Eq. (2).

Further insights and numerical solutions for Eqs. (1), (2) and (3) can be found in Basawa et al. (2004) and Lit (2020).

Predictions were generated using the Time Series Lab (Lit et al. 2022), Score Edition software V. 1.5 (https://timeserieslab.com), assisted by Visual Recurrence Analysis (https://visual-recurrence-analysis.software.informer.com/4.9) and CurveExpert Professional 1.6 (https://www.curveexpert.net).
## 2.5 Exogenous input model structure

In this study, we used exogenous climatic input to predict and forecast NSD time-series. Considering the known natural temporal variabilities of NSD, devoid of anthropogenic influence as noted by Diodato et al. (2023), the PARMAX(TVAR)-EGARCH model was selected. This model was chosen due to its ability to represent the natural variability in NSD and effectively respond to exogenous climate forcing.

To capture the combined influence of the atmospheric pressure and temperature on climate patterns and support the input data model, we introduced a dedicated dimensionless index, the Teleconnection-Climate Pattern Index (TCPI), was introduced. The TCPI is defined as:

$$TCPI = \frac{1}{(AO + Tw + 6)} \qquad (4)$$

This unique, dimensionless index is designed to account for both teleconnection patterns, (through the Arctic Oscillation, AO, unitless) and local climate conditions (through the mean winter temperature, Tw, °C). The constant “6” is an empirical adjustment term included to balance the influence of Tw within the TCPI equation. Notably, the variable X in Eq. (2) must encompass the TCPI.

In the TCPI, the AO and Tw values are not explicitly weighted beyond the standardisation of AO. The AO, in line with the findings of Capozzi et al. (2022), is used as a teleconnection index of the large-scale atmospheric patterns. On the other hand, Tw serves as an indicator of small-scale thermal forcing, as highlighted by Lehr et al. (2012). This formulation thus ensures that the TCPI provides a robust metric for investigating the combined effects of large-scale atmospheric variability and localised temperature variations on the climate system.

The inclusion of TCPI features in the PARMAX(TVAR)-EGARCH model helped to capture the behaviour of an autoregressive time-series data, reflecting the dynamic interplay between atmospheric pressure and temperature variations. Importantly, the ratio of the TCPI term in Eq. (4) reveals a stronger relationship with the NSD when the exogenous variables
AO and Tw are examined individually. Here, Tw denotes the mean winter temperature at Montevergine from December to April. This approach ensures a comprehensive understanding of the climate pattern by considering both large- and small-scale indices. It should be considered as a straightforward support indicator with two variables incorporated in the model, which gave satisfactory results while limiting the uncertainties in the information flow of the PARMA(TVAR)-EGARCH model, ensuring robustness and reliability.

Exogenous input variables were obtained from different monthly sources. The Arctic Oscillation (AO) index in our study was originally obtained from Trenberth and Paolino (1980) and subsequently updated using Climate Explorer (https://climexp.knmi.nl). The AO index is a dimensionless, standardised anomaly computed from the Geopotential Height (GPH) data at the 500 and 850 hPa pressure levels. It is derived by calculating the zonal mean area-weighted difference between the 30° N-60° N latitude band and the 60° N-90° N latitude band. This difference captures the variability in atmospheric pressure patterns between the mid-latitudes and the Arctic, influencing general circulation and climate patterns over the Northern Hemisphere. The AO index is expressed in standardised units (z-scores), representing the deviation from its long-term mean, divided by the standard deviation of the index over a reference period. This normalisation ensures that the AO is dimensionless and allows for consistent comparisons across different time periods and studies. Positive AO values correspond to lower-than-mean pressure over the Arctic and higher-than-mean pressure over the mid-latitudes, while negative AO values indicate the opposite.

The mean winter temperature (Tw) used in our study reflects the mean temperature from December to April and is expressed in degrees Celsius (°C). The Tw data for Montevergine were extracted from the original registers of the Montevergine Observatory, as documented
by Carbone (2014). These values are used directly without normalisation, preserving their
original units (°C). Data updates until 2022 were performed using the E-OBS temperature
dataset from the Copernicus Climate Change Service (C3S,
https://surfobs.climate.copernicus.eu), accessed via Climate Explorer (Cornes et al. 2018).

**3. Results and Discussion**

Temporal prediction features operate across multiple time scales, reflecting the inherent
dual nature of snowfall day processes that involve both short- and long-term memory.
However, the existence of multiple scales poses a substantial challenge in detecting,
interpreting and conveying different aspects of snowfall day systems, as well as managing
their forecasts. Exogenous variables play a crucial role in providing supportive information
in this complex process, ultimately enhancing the effectiveness of the forecasting phase.

**3.1 Influence of the exogenous (climate) variables**

In this study, we employed the PARMAX(TVAR)-EGARCH model to investigate and
analyse the direct influence of an exogenous variable, considering a paired set of variables to
account for the influence of external climatic forcings in decadal projections. Specifically, we
used exogenous data from the Teleconnection-Climate Pattern Index (TCPI), which combines
large-scale factors like the Arctic Oscillation with smaller-scale factors such as air
temperature, depicting major climate variability in the central Mediterranean area.

Figure 7a depicts the linear relationship between the TCPI index and the annual NSD. A
notable observation is that only a few data points fall within the 95% prediction boundaries,
with just eight points falling outside. The p-value in the ANOVA is less than 0.05, with
Pearson's correlation coefficient equal to 0.74, indicating a statistically significant
relationship between the number of snowfall days and the teleconnection climate pattern index at the 95% confidence level.

To assess the individual contributions of both the AO and Tw indicators to the TCPI index, we examined the highest *p*-value on the independent variables, which is 0.00358 and corresponds to AO. Since this *p*-value is less than 0.05, the AO term is statistically significant at the 95% confidence level. Consequently, there is an indication that we should retain all variables in the TCPI information. Subsequently, the TCPI index time-series for the period 2023-2060 was explored (see Appendix), and this information was then incorporated into the flow of NSD information.

## 3.2 Competitor model

The selection of the best-performing model for the NSD time-series involved the calibration of autoregressive parameters. Leveraging the suitability of the ARMA-EGARCH Time Series Lab model family for hydrological and periodic time-series, these models were employed to estimate the parameters of the competing model. Ultimately, the PARMAX(TVAR)-EGARCH model was chosen for the NSD time-series due to its proficiency in recognising non-stationary signals and effectively incorporating information from exogenous factors. The autoregressive parameters in the PARMAX(TVAR)-EGARCH model were carefully selected to minimise the model's RMSE, MAE and MAPE. In this model with a Gaussian distribution function, the autoregressive component is non-stationary and the sum of the *m* periodic components must equal zero. Achieving this involves estimating the first *m*−1 components, determining the final component through a combination with the zero-sum limit (Lit 2020). The frequency spectrum graphically shows the length of *m* and the amplitude of cyclic components within a specific frequency range.
The frequency spectrum for the PARMAX(TVAR)-EGARCH periodic model parameters
fitted to the NSD data (detrended and standardised) is given in Figure 7b. Notably, the main
cycles are observed between 1/f(0.01) and 1/f(0.1) (Figure 7b, red box). Accordingly, the
cycle was set to 1/f(0.040), corresponding to a period *m* of about 28 years, strategically
chosen within the high cycle window and supported by validation (after several trials). This
aligns with the observed trend in the time-series plot during both the training and validation
periods (Figure 7c, light blue and orange, respectively).

3.3 *Model parameterisation and assessment*

As indicated in Figure 7c, the development of observed and modelled data is presented using
the PARMAX-EGARCH model across the training period (1884-1991) and validation period
(1992-2022). Notably, the timeline in the validation period aligns closely with that of the
training epoch, evident in the agreement between the observation (blue) and model (orange)
lines. The error statistics for both the training and validation periods are summarised in Table
1. Model validation error statistics were computed using the last 30 data values in the time-
series, emphasising the need to confirm forecast accuracy for reliable predictions. The
consistency in outcomes between the training period (1884-1991, Table 1) and the validation
period (1992-2022) suggests that the model is likely to perform as predicted in future
forecasts.

Furthermore, a comparison was conducted between the PARMAX(TVAR)-EGARCH and
PARMA models during validation, with PARMA lacking consideration of exogenous support
variables. A scrutiny of the second and third columns in Table 1 reveals that the accepted
model, incorporating external input, was an appropriate strategy for predicting NSD.

There is no established cutoff point for determining what constitutes an appropriate Log Loss
value. The Log Loss function, with values between 7 and 8, serves as a measure of uncertainty.
A MAPE of about 20%, on the other hand, stands out, indicating a reasonable degree of accuracy (<25%, Swanson et al. 2011) for projections aiming to capture both the long-term trend and pattern of NSD. Additionally, the Pearson's correlation coefficient and Kling-Gupta Efficiency index, both substantially exceeding 0.5, indicate that the model outperforms the benchmark predictor represented by the mean of observations, with minimal uncertainty in the model's estimations. The parameterisation of the PARMAX(TVAR)-EGARCH model for forecasting was conducted throughout the whole NSD time-series from 1884 to 2022 (139 years), ensuring that the most recent portion of the time-series is available as valuable information for future predictions.

Given the PARMAX [1,1] [1,1]28 - EGARCH[1,1]3 relationship in the NSD model, as proven by the training predictions, the model is effectively implemented, and its values are optimised by refining process parameters for convergence via numerical derivatives (log-likelihood=-466). The model's appropriateness was tested at a significance level of 0.05 using the constant residuals assumption.

## 3.4 Normality diagnostic checks

We can now assess the normality of the PARMAX(TVAR)-EGARCH model through various goodness-of-fit tests to determine if the residuals form a white noise pathway (Jere and Moyo 2016). The histogram of Figure 8a shows that the residual distribution is Gaussian, with the mean residual approaching zero. Normality can be assessed using QQ plots. In Figure 8b, the QQ plot of the residuals indicates a normal distribution as the points (empty black dots) closely align with the blue straight line within the confidence interval (blue dashed line).

## 3.5 Number of snowfall days forecast and large-scale climate variability impact
Figure 9 presents a visual analysis of annual projections for NSD (orange line), revealing that the extremes of future value oscillations are likely to remain within the range of variability (±1σ) observed over the last 30 years (blue and orange thick lines). Notably, important interdecadal oscillations are evident. Consequently, the conclusion of the current decade is anticipated to continue experiencing NSD levels above the mean value (thin gray line). Subsequent forecasts suggest a decline in expectations during the 2030s and 2040s, with NSD approaching -1σ between 2034 and 2038, and between 2045 and 2049. Especially in the final period, a positive shift is observed, starting from 2050, indicating a mean increase from about 20 to 30 snowfall days per year. During these times, spring seasons would witness a change toward pluvial-dominated sediment action due to less snow and greater erosive rainfall (Zhang et al. 2023). Periods with NSD exceeding the mean values are expected within the years 2039-2041 and 2051-2057. The observed oscillations, particularly beyond 2030, are consistent with findings from the Euro-CORDEX models (Faggian et al. 2024) under the emission pathways RCP8.5 (without mitigation). These models project that the Alto Adige region is expected to encounter increased wet snow loads (WSL) in the winter season over the next 40 years (Figure 9, dashed gray line). This is noteworthy, especially in contrast to the projected reductions in WSL, at lower sites in Abruzzo (central Italy) and Frabosa (northwestern Italy).

The forecasted interdecadal oscillations in NSD at Montevergine are accompanied by overlapping interannual variations, although these variations are notably less pronounced than during the historical period. The long-term trend suggests a subtle increase starting around the year 2030. The results point out a relatively complex pattern of climate change effects: the number of snowfall days, considered in both their interannual and interdecadal variations, as well as their long-term trend, is expected to generally increase. However, it remains unclear whether this increase will affect all elevations or primarily higher ones, given the potential
shift from snowfall to rainfall if global warming persists. Nevertheless, considering Montevergine’s location in the southern Apennines, an intensification of NSD is not only anticipated over the higher Apennine regions but also aligns with the forecast for Alto Adige on both interdecadal and longer-term scales.

Following recommendations to convey information on the uncertainty of scenario forecasts, we examined the agreement among the PARMAX(TVAR)-EGARCH models with respect to weather patterns across the European region. To achieve this, we considered a set of 22 ensemble simulations from Euro-CORDEX: 12 high-resolution (~12 km) regional models (https://www.euro-cordex.net) with the RCP8.5 emissions scenario, from which a climatic wet snow accretion model named CLIWAC was derived (Faggian et al. 2024).

The Sixth Coupled Model Intercomparison Project (CMIP6) projects an increase in percentiles and the potential for expected extreme magnitudes of daily snowfall across global land regions until at least the middle of the 21st century (Quante et al. 2021). However, Croce et al. (2018) showed that these increases do not necessarily result in a decrease of snow loads. Regarding the impact of large-scale climate variability impact on snowfall, Jeong et al. (2022) identified a significant positive relationship between tropical convection activity and the AO in the reanalysis, suggesting that increased sea-surface temperature (SST)-forced convection over the Indian Ocean (IO) could lead to a positive trend in the AO. Consequently, the AO is projected to exhibit a growing trend in the future climate, given the anticipated IO warming phase in the coming decades. However, there is substantial uncertainty about future changes in IO frequency and intensity (Endris et al. 2019), exacerbated by challenges in modelling dynamic changes in atmospheric circulation (Smith et al. 2020).

Climate models generally show satisfactory agreement on long-term temperature signals (Smith et al. 2020), with a decadal mean global temperature range of 1.3-3.4 °C projected by 2050 (Lehner et al. 2023). Nevertheless, some studies (e.g. Otto et al. 2013; Stott et al. 2013)
suggest that models with the highest transient climate response may overestimate near-term
warming (here, 2016-2035). This highlights an area of ongoing research in refining climate
models for improved accuracy across different timescales. In fact, this discrepancy persists in
CMIP6 models, showing less consistency with observed temperature trends over previous
decades (Bonnet et al. 2021). An additional consideration is the potential impact of Arctic
amplification on Northern Hemisphere terrestrial snow, indicating potential increases in the
Arctic/Siberian snowpack coinciding with sea-ice loss (Ghatak et al. 2010). These effects may
intensify in future decades with ongoing global warming (Henderson et al. 2018), potenitally
enhancing the linkage between the snow-AO teleconnection and creating more favourable
conditions for snowfall in southern Europe. However, a recent study by Cohen et al. (2023)
suggests that the impact of Arctic changes on cold extreme weather in southern Europe is
complex and poorly understood, highlighting the need for more in-depth investigation into the
physical origins of these observed connections. Furthermore, debates persist regarding the
non-stationary or non-causal nature of snow teleconnection, fuelling ongoing discussions
around their long-term reliability (Henderson et al. 2018; Peings et al. 2013).

3.6 Cautions and errors in forecasting assessment

Even when an accurate model has been generated with both low-bias and low variance the forecasts errors are likely to increase with the forecast horizon, given there is usually an interdependency between values which are close in time (Haben et al. 2023). Hence, if comparing the errors for forecasting tomorrow, versus forecasts for the following day and so on, there would be an expected upward trajectory in the forecast errors/scores, but it may not be that simple. This is because there is more variation in demand (and hence larger errors on average) during certain periods of the day compared to others. However, even in this case, average daily errors appear to be increasing.
Another source of forecast error may be due to their dependence on the input variables. Many load forecasts are strongly related to weather and therefore weather forecasts are utilised within the load forecast models. However, if these inputs are inaccurate (for example through measurement, forecasts or even calibration errors) then the load forecast will also be inaccurate. In other words, for surprising or unusual errors in the data, the input variables should also be considered as possible causes. In addition, external forcing acting on the model introduces additional sources of error, in particular from the divergence between the increase in mean temperature and snowfall (Quante et al. 2021). The rise in temperature reduces the likelihood of snowfall in Western and Southern Europe, in line with both past observations (Fischer and Knutti, 2016) and future projections (Prein and Heymsfield 2020). The expected enhancement of the hydrological cycle is predicted to transtate into more rain than snowfall, further contributing to forecast inaccuracies (Quante et al. 2021). It is important to recognise and address these complexities in the assessment of projections in order to improve the robustness and reliability of models.

# 4. Conclusions and future directions

## 4.1 Novel statistical approaches in decadal climate predictions

In the dynamic field of decadal prediction, which increasingly relies on models of intermediate complexity, our study addresses the challenges associated with statistical modeling, particularly regarding precipitation patterns in Mediterranean regions in response to climate change. The limited progress in this area highlights the importance of fostering communication and collaboration between scientists. Our research seeks to bridge this gap by introducing a quite novel statistical forecasting approach that links observed snowfall patterns (annual time-series) with historical changes and extrapolates these patterns into the future.
Specifically, our study proposes a data-driven, seasonal (periodic) autoregressive statistical technique to generate 38-year (2022-2060) estimates of snowfall days for the Montevergine site. Based on local historical observations from 1884 to 2021, our approach provides a robust and nuanced forecasting method that improves our understanding of multidecadal climate variability trends. This method offers valuable insights for climate science research and practical applications. As we navigate the ever-evolving landscape of climate studies, the integration of statistical techniques, as exemplified by the approach presented in this study, is emerging as a critical step toward achieving more accurate and informed decadal predictions. By helping to refine prediction methods, our work contributes to the advancement of climate science and provides researchers, policymakers and stakeholders with improved tools to understand and prepare for the complex challenges posed by climate change.

Additionally, this study explored the intricate landscape of decadal snowfall forecasting, providing nuanced insights into the complex interplay between climate variables, large-scale patterns and temporal forecasting challenges. The use of the PARMAX(TVAR)-EGARCH model stands out as a notable contribution, demonstrating its effectiveness in capturing the dual nature of snowfall processes, incorporating both short- and long-term memory. The robustness of the chosen model, evidenced by its consistent performance throughout the training and validation periods, underlines a notable strength.

## 4.2 Incorporating teleconnection indices

Chase et al. (2006) provide insights into the importance and utility of teleconnection indices in climate research. In our study, we use teleconnection indices, which are pivotal in integrating diverse forcing agents and providing a comprehensive perspective on the influences of various factors on climatic changes. Teleconnection indices effectively capture the combined effects of multiple forcing agents, including anthropogenic influences, on
climate patterns. It is important to note that our predictions are valid under the current forcing scenario. In our statistical model forecast, the direct correlation of an observational time-series of global greenhouse gases or aerosols with a site-scale precipitation record is unlikely to improve predictive power. Instead, this integrative approach has been shown to yield robust responses beyond insights gained from considering individual forcing agents in isolation. If emissions were to change significantly in the coming years, it would be necessary to update our forecasts accordingly. However, this requirement for updates is not a specific weakness of our approach. The Intergovernmental Panel on Climate Change (IPCC) regularly provides assessments to update the understanding of future climate trends. Continuous monitoring and model adjustments are thus essential for gaining a better understanding of the direction our climate is heading. Even within the context of time-series analysis, predictions can be repeated using the same methodology. Detecting departures from previously predicted patterns can provide insights into possible local-scale changes that climate models may not be able to capture.

## 4.3 Physical significance of periodicities

Periodic behaviour in climatic time-series is a well-documented phenomenon. Periodicities in climatic data often reflect underlying physical processes (e.g., Burroughs 2003), which can be influenced by a combination of natural cycles, atmospheric oscillations and other environmental factors. For example, natural cycles such as the 11-year sunspot cycle, atmospheric oscillations like the Arctic Oscillation (AO) and the North Atlantic Oscillation (NAO), and periodic changes in sea surface temperatures associated with the Atlantic Multidecadal Oscillation (AMO) can impose periodic behaviour on climate-related data. The IPCC AR5 WG1 report (Stocker et al., 2013) discusses these influences on climate: the AO and the NAO (Chapter 2), sunspot cycle variability (Chapter 8) and the AMO (Chapter 14).
Each geographic location can exhibit unique periodicities due to its specific climatic and environmental influences. The identification of a 28-year period within the NSD (number of snowfall days) time-series data from Montevergine has significant implications for understanding the region's climatic patterns. However, in the context of snow cover variability, different periodicities have been identified in other regions. For instance, Valt and Cianfarra (2010) reported several dominant periodicities in the Alps, namely 39.5, 19.7, 11.2, and 7.2 years. These periods correspond to different climatic drivers and oscillatory processes affecting that particular region. While the identified periods differ from the 28-year cycle found in Montevergine, they illustrate the presence of multiple periodicities in snowfall data driven by regional climatic dynamics.

### 4.4 Integration of broader climate dynamics and future research

The integration of exogenous variables, in particular the Teleconnection-Climate Pattern Index (TCPI), enriches the prediction capabilities and highlights the importance of considering broader climate dynamics in snowfall forecasting. Numerous tests were carried out to include various exogenous variables, but only the inclusion of TCPI yielded consistently positive results. The introduction of variables, such as AO and Tw, gave satisfactory results and improved the robustness of the model. However, minimal improvements from additional variables are often redundant and may introduce further uncertainties into the forecast horizon. These uncertainties need to be managed within the information flow of the main PARMAX(TVAR)-EGARCH model.

Attention to normality diagnostics and careful calibration of model parameters lend confidence to the results of the study. However, given the inherent complexity of climate systems and uncertainties in long-term projections, the focus on the Montevergine Observatory, while providing valuable regional insights, may not fully capture the diversity of
climatic conditions in other geographical settings. With these limitations in mind, looking
ahead, there is speculation about the observed changes in snow days, with an expected decline
in the 2030s and 2040s, raising questions about wider implications for ecosystems and water
resource management.

Future research priorities should include exploring the intricacies of the effects of climate
variability on snowfall, particularly the observed links with large-scale patterns such as the
Arctic Oscillation. Broadening the geographical scope of the study to include different
terrains and incorporating advances in climate modelling are important ways of improving the
generalisability of the results and addressing uncertainties in global warming projections.

The complexity of snow teleconnections and ongoing debates about their non-stationary or
non-causal nature emphasise the need for a comprehensive understanding. This study
contributes to the ongoing discourse by shedding light on specific connections while
highlighting the intricate, multifaceted nature of these relationships. The evolving field of
Arctic change and its implications for southern Europe requires continued scrutiny and
collaboration among researchers, highlighting the importance of an open and collaborative
approach to advancing our understanding of climate dynamics.

In summary, while our approach has limitations and acknowledges the need for ongoing
monitoring and updates, incorporating teleconnection indices and continuously evaluating
predictions can enhance our understanding of future climate changes beyond what individual
forcing agents can provide.

Declarations

Funding

This was an investigator-driven study with no financial support.

**Disclosure statement:** No potential conflict of interest was reported by the authors.
**Authors' contributions:** Conceptualisation, methodology design and preparation of figures were carried out by N.D. S.R. put the modelling approach into the context of regional hydrology. G.B. supervised the study. All co-authors reviewed and contributed to the final version of the manuscript.

<figure><img src="image_11.png" /></figure>
References

Azari, A. et al., 2021. Integrated preprocessing techniques with linear stochastic approaches in groundwater level forecasting. *Acta Geophysica*, 69, 1395-1411.

https://doi.org/10.1007/s11600-021-00617-2

Bailey, H. et al., 2021. Arctic sea-ice loss fuels extreme European snowfall. *Nature Geoscience*, 14, 283–288. https://doi.org/10.1038/s41561-021-00719-y

Baltaci, H. et al., 2020. Long-term variability and trends of extended winter snowfall in Turkey and the role of teleconnection patterns. *Meteorological Applications* 27, e1891.

https://doi.org/10.1002/met.1891

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

https://doi.org/10.1038/nature04141

Basawa, I.V., Lund, R., Shao, Q., 2004. First-order seasonal autoregressive processes with periodically varying parameters. *Statistics & Probability Letters*, 67, 299–306.

https://doi.org/10.1016/j.spl.2004.02.001

Bednorz, E., 2004. Snow cover in eastern Europe in relation to temperature, precipitation and circulation. *International Journal of Climatology*, 24, 591–601.

https://doi.org/10.1002/joc.1014

Beniston, M., 1997. Variations of snow depth and duration in the Swiss Alps over the last 50 years: Links to changes in large-scale climatic forcings. *Climatic Change*, 36, 281–300.

https://doi.org/10.1023/A:1005310214361

Beniston, M. et al., 2014. Assessing climate change impacts on the quantity of water in Alpine regions: Foreword to the adaptation and policy implications of the EU/FP7 "ACQWA" project. *Environmental Science & Policy*, 43, 1–4.

https://doi.org/10.1016/j.envsci.2014.01.009
Beniston, M. et al., 2018. The European mountain cryosphere: a review of its current state, trends, and future challenges. *The Cryosphere*, 12, 759–794. https://doi.org/10.5194/tc-12-759-2018

Blahušiaková, A. et al., 2020. Snow and climate trends and their impact on seasonal runoff and hydrological drought types in selected mountain catchments in Central Europe. *Hydrological Sciences Journal*, 65, 2083–2096.
https://doi.org/10.1080/02626667.2020.1784900

Bonnet, R. et al., 2021. Increased risk of near term global warming due to a recent AMOC weakening. *Nature Communications*, 12, 6108. https://doi.org/10.1038/s41467-021-26370-0

Box, G.E.P., Jenkins, G.M., Reinsel, G.C., 2013. *Time series analysis: Forecasting and control*. Hoboken, NJ: Wiley.

Brázdil, R. et al., 2023. Severity of winters in the Czech Republic during the 1961–2021 period and related environmental impacts and responses. *International Journal of Climatology*, 43, 2820–2842. https://doi.org/10.1002/joc.8003

Burroughs, W.J., 2003. *Weather cycles: real or imaginary?* Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511841088

Capozzi, V., Budillon, G., 2017. Detection of heat and cold waves in Montevergine time series (1884–2015). *Advances in Geosciences*, 44, 35–51. https://doi.org/10.5194/adgeo-44-35-2017

Capozzi, V., De Vivo, C., Budillon, G., 2022. Synoptic control over winter snowfall variability observed in a remote site of Apennine Mountains (Italy), 1884–2015. *The Cryosphere*, 16, 1741–1763. https://doi.org/10.5194/tc-16-1741-2022

Carbone, M., 2014. *Ricostruzione della serie storica delle temperature registrate all'Osservatorio di Montevergine*. Thesis. University of Sannio. (in Italian)
Chase, T.N., Pielke Sr, R.A., Avissar, R., 2006. Teleconnections in the Earth System. In: Anderson, M.G., McDonnell, J.J., eds. *Encyclopedia of Hydrological Sciences*, 15. *Global Hydrology*. Hoboken, NJ: John Wiley & Sons, 183.
https://doi.org/10.1002/0470848944.hsa190

Cohen, J. et al., 2023. No detectable trend in mid-latitude cold extremes during the recent period of Arctic amplification. *Communications Earth & Environment*, 4, 1–9.
https://doi.org/10.1038/s43247-023-01008-9

Cornes, R.C. et al., 2018. An ensemble version of the E-OBS temperature and precipitation data sets. *Journal of Geophysical Research*, 123, 9391–9409.
https://doi.org/10.1029/2017JD028200

Croce, P. et al., 2018. The snow load in Europe and the climate change. *Climate Risk Management*, 20, 138–154. https://doi.org/10.1016/j.crm.2018.03.001

Diodato, N., Ljungqvist, F.C., Bellocchi, G., 2021. Empirical modelling of snow cover duration patterns in complex terrains of Italy. *Theoretical and Applied Climatology*, 147, 1195–1212. https://doi.org/10.1007/s00704-021-03867-8

Diodato, N. et al., 2023. Climatic fingerprint of spring discharge depletion in the southern Italian Apennines from 1601 to 2020 CE. *Environmental Research Communications*, 4, 125011. https://doi.org/10.1088/2515-7620/acae23

Diodato, N., Büntgen, U., Bellocchi, G., 2019. Mediterranean winter snowfall variability over the past millennium. *International Journal of Climatology*, 39, 384–394.
https://doi.org/10.1002/joc.5814

Diodato, N., Fratianni, S., Bellocchi, G., 2020. Reconstruction of snow days based on monthly climate indicators in the Swiss pre-alpine region. *Regional Environmenta Change*, 20, 55. https://doi.org/10.1007/s10113-020-01639-0
Diodato, N., Lanfredi, M., Bellocchi, G., 2023. Long-range, time-varying statistical prediction of annual precipitation in a Mediterranean remote site. *Environ. Res.: Climate* 2, 045011.
https://doi.org/10.1088/2752-5295/acffe9

Dye, D.G., 2002. Variability and trends in the annual snow-cover cycle in Northern Hemisphere land areas, 1972–2000. *Hydrological Processes*, 16, 3065–3077.
https://doi.org/10.1002/hyp.1089

Efthymiadis, D. et al., 2007. Influence of large-scale atmospheric circulation on climate variability in the Greater Alpine Region of Europe. *Journal of Geophysical Research*, 112, D12104. https://doi.org/10.1029/2006JD008021

Endris, H.S. et al., 2019. Future changes in rainfall associated with ENSO, IOD and changes in the mean state over Eastern Africa. *Climate Dynamics*, 52, 2029–2053.
https://doi.org/10.1007/s00382-018-4239-7

Enzi, S., Bertolin, C., Diodato, N., 2014. Snowfall time-series reconstruction in Italy in the last 300 years. *The Holocene*, 24, 346-356. http://hol.sagepub.com/content/24/3/346.short

Faggian, P., Trevisiol, A., Decimi, G., 2024. Future projections of wet snow frequency and wet snow load on overhead high voltage conductors over Italy. *Cold Regions Science and Technology*, 217, 103980. https://doi.org/10.1016/j.coldregions.2023.103980

Fischer, E.M., Knutti, R., 2016. Observed heavy precipitation increase confirms theory and early models. *Nature Climate Change*, 6, 986–991. https://doi.org/10.1038/nclimate3110

Fratianni, S. et al., 2015. How snow and its physical properties change in a changing climate alpine context? In: Lollino, G. et al., eds. *Engineering Geology for Society and Territory*. Cham: Springer International Publishing, 1, 57-60. https://doi.org/10.1007/978-3-319-09300-0_11

Fuente, I. de la, Navarro, E., Serna, G., 2020. Reverse Mortgage Risks. Time Evolution of VaR in Lump-Sum Solutions. *Mathematics* 8, 2043. https://doi.org/10.3390/math8112043
Gazzolo, T., 1973. Le precipitazioni nevose in Italia. *Bollettino della Società Geografica Italiana*, 2, 11-33. Retrieved from [https://bsgi.it/index.php/bsgi/article/view/6145](https://bsgi.it/index.php/bsgi/article/view/6145). (in Italian)

Ghatak, D. et al., 2010. On the emergence of an Arctic amplification signal in terrestrial Arctic snow extent. *Journal of Geophysical Research*, 115, D24105.

[https://doi.org/10.1029/2010JD014007](https://doi.org/10.1029/2010JD014007)

Haben, S., Voss, M., Holderbaum, W., 2023. Verification and evaluation of load forecast models. In: Haben, S., Voss, M., Holderbaum, W., eds. *Core concepts and methods in load forecasting: with applications in distribution networks*. Cham, Springer International Publishing, 89–105. [https://doi.org/10.1007/978-3-031-27852-5_7](https://doi.org/10.1007/978-3-031-27852-5_7)

Hale, K.E. et al., 2023. Recent decreases in snow water storage in western North America. *Communications Earth & Environment*, 4, 1–11. [https://doi.org/10.1038/s43247-023-00751-3](https://doi.org/10.1038/s43247-023-00751-3)

Hamouda, M.E., Pasquero, C., Tziperman, E., 2021. Decoupling of the Arctic Oscillation and North Atlantic Oscillation in a warmer climate. *Nature Climate Change*, 11, 137–142.

[https://doi.org/10.1038/s41558-020-00966-8](https://doi.org/10.1038/s41558-020-00966-8)

Henderson, G.R. et al., 2018. Snow-atmosphere coupling in the Northern Hemisphere. *Nature Climate Change*, 8, 954–963. [https://doi.org/10.1038/s41558-018-0295-6](https://doi.org/10.1038/s41558-018-0295-6)

Jere, S., Moyo, E., 2016. Modelling epidemiological data using Box-Jenkins procedure. *Open Journal of Statistics*, 6, 295–302. [https://doi.org/10.4236/ojs.2016.62025](https://doi.org/10.4236/ojs.2016.62025)

Jeong, Y.-C. et al., 2022. Indian Ocean warming as key driver of long-term positive trend of Arctic Oscillation. *npj Climate and Atmospheric Science*, 5, 1–10.

[https://doi.org/10.1038/s41612-022-00279-x](https://doi.org/10.1038/s41612-022-00279-x)
Kis, A., Pongrácz, R., 2021. The role of temperature and the NAO index in the changing snow-related variables in European regions in the period 1900–2010. *Hungarian Geographical Bulletin*, 70, 325–337. https://doi.org/10.15201/hungeobull.70.4.3

Klein, G. et al., 2016. Shorter snow cover duration since 1970 in the Swiss Alps due to earlier snowmelt more than to later snow onset. *Climatic Change*, 139, 637–649.
https://doi.org/10.1007/s10584-016-1806-y

Koopman, S.J., Lucas, A., Scharth, M., 2016. Predicting time-varying parameters with parameter-driven and observation-driven models. *The Review of Economics and Statistics*, 98, 97–110. https://doi.org/10.1162/REST_a_00533

Lehr, C., Ward, P.J., Kumma, M., 2012. Mountain Research and Development 32, 431-445.

Lehner, F. et al., 2023. New potential to reduce uncertainty in regional climate projections by combining physical and socio-economic constraints. *AGU Advances*, 4, e2023AV000887.
https://doi.org/10.1029/2023AV000887

Lit, R., 2020. *Forecasting the VIX in the midst of COVID-19* [online]. Available from:
https://timeserieslab.com/articles/rlit_vix_v2.pdf

Lit, R., Koopman, S.J., Harvey, A.C., 2022. *Time Series Lab - Advanced Time Series Forecasting Software* [online]. Available from: https://timeserieslab.com

Mandelbrot, B.B., Taqqu, M.S., 1979. Robust R/S analysis of long run serial correlation.
*Bulletin of the International Statistical Institute*, 48, 59-104.

Marcolini, G. et al., 2017. Variability in snow depth time series in the Adige catchment.
*Journal of Hydrology: Regional Studies*, 13, 240-254.
https://doi.org/10.1016/j.ejrh.2017.08.007

Martines, N.M.L., 2018. *La memoria: un tema letterario attraverso i secoli* [online].
Available from:
https://www.treccani.it/magazine/chiasmo/lettere_e_arti/Memoria/memoria_ssc_infinito_mondo_della_memoria.html

Marty, C., Tilg, A.-M., Jonas, T., 2017. Recent Evidence of Large-Scale Receding Snow Water Equivalents in the European Alps. *Journal of Hydrometeorology*, 18, 1021–1031.
https://doi.org/10.1175/JHM-D-16-0188.1

McLeod, A.I., Li, W.K., 1983. Diagnostic Checking Arma Time Series Models Using Squared-Residual Autocorrelations. *Journal of Time Series Analysis*, 4, 269–273.
https://doi.org/10.1111/j.1467-9892.1983.tb00373.x

Meredith, M. et al., 2019. Polar regions. In: Pörtner, D.C. et al., eds. *IPCC Special Report on the Ocean and Cryosphere in a Changing Climate*. Cambridge and New York, NY: Cambridge University Press, 3, 203–320. https://doi.org/10.1017/9781009157964.005

Merino, A. et al., 2014. Snowfall in the Northwest Iberian Peninsula: Synoptic circulation patterns and their influence on snow day trends. *The Scientific World Journal*, 2014, e480275. https://doi.org/10.1155/2014/480275

Misra, A. et al., 2020. Topographic and climatic influence on seasonal snow cover: Implications for the hydrology of ungauged Himalayan basins, India. *Journal of Hydrology*, 585, 124716. https://doi.org/10.1016/j.jhydrol.2020.124716

Montgomery, D.C., Jennings, C.L., Kulahci, M., 2015. Introduction to time series analysis and forecasting. Hoboken, NJ: John Wiley & Sons.

Musselman, K.N. et al., 2021. Winter melt trends portend widespread declines in snow water resources. *Nature Climate Change*, 11, 418–424. https://doi.org/10.1038/s41558-021-01014-9

Nerantzaki, S.D., Hristopulos, D.T., Nikolaidis, N.P., 2020. Estimation of the uncertainty of hydrologic predictions in a karstic Mediterranean watershed. *Science of The Total Environment*, 717, 137131. https://doi.org/10.1016/j.scitotenv.2020.137131
Nigrelli, G., Collimedaglia, M., 2012. Reconstruction and analysis of two long-term
precipitation time series: Alpe Devero and Domodossola (Italian Western Alps).
*Theoretical and Applied Climatology*, 109, 397-405. https://doi.org/10.1007/s00704-012-0586-1

Notarnicola, C., 2020. Observing snow cover and water resource changes in the high
mountain Asia region in comparison with global mountain trends over 2000–2018. *Remote
Sensing*, 12, 3913. https://doi.org/10.3390/rs12233913

Ostertagova, E., Ostertag, O., 2012. Forecasting using simple exponential smoothing method.
*Acta Electrotechnica et Informatica*, 12, 62–66. https://doi.org/10.2478/v10198-012-0034-2

Otto, A. et al., 2013. Energy budget constraints on climate response. *Nature Geoscience*, 6, 415–416. https://doi.org/10.1038/ngeo1836

Peings, Y. et al., 2013. How stationary is the relationship between Siberian snow and Arctic Oscillation over the 20th century? *Geophysical Research Letters*, 40, 183–188.
https://doi.org/10.1029/2012GL054083

Pinna, M., 1973. La durata del manto nevoso in Italia. *Bollettino della Società Geografica Italiana*, 2, 35-62. Retrieved from https://bsgi.it/index.php/bsgi/article/view/6146 (in Italian)

Prein, A.F., Heymsfield, A.J., 2020. Increased melting level height impacts surface
precipitation phase and intensity. *Nature Climate Change*, 10, 771–776.
https://doi.org/10.1038/s41558-020-0825-x

Quante, L. et al., 2021. Regions of intensification of extreme snowfall under future warming.
*Scientific Reports*, 11, 16621. https://doi.org/10.1038/s41598-021-95979-4
Rubinetti, S. et al., 2020. Robust decadal hydroclimate predictions for Northern Italy based on a twofold statistical approach. *Atmosphere*, 11, 671.

https://doi.org/10.3390/atmos11060671

Sandgathe, S. et al., 2020. Exploring the need for reliable decadal prediction. *Bulletin of the American Meteorological Society*, 101, E141–E145. https://doi.org/10.1175/BAMS-D-19-0248.1

Sarhadi, A., Kelly, R., Modarres, R., 2014. Snow water equivalent time-series forecasting in Ontario, Canada, in link to large atmospheric circulations. *Hydrological Processes* 28, 4640–4653. https://doi.org/10.1002/hyp.10184

Schroeer, K., Kirchengast, G., 2018. Sensitivity of extreme precipitation to temperature: the variability of scaling factors from a regional to local perspective. *Climate Dynamics*, 50, 3981–3994. https://doi.org/10.1007/s00382-017-3857-9

Shi, S., Liu, G., 2021. The latitudinal dependence in the trend of snow event to precipitation event ratio. *Scientific Reports*, 11, 18112. https://doi.org/10.1038/s41598-021-97451-9

Smith, D.M. et al., 2020. North Atlantic climate far more predictable than models imply. *Nature*, 583, 796–800. https://doi.org/10.1038/s41586-020-2525-0

Smith, T., Bookhagen, B., 2018. Changes in seasonal snow water equivalent distribution in High Mountain Asia (1987 to 2009). *Science Advances*, 4, e1701550.

https://doi.org/10.1126/sciadv.1701550

Solaraju-Murali, B. et al., 2021. Multi-annual prediction of drought and heat stress to support decision making in the wheat sector. *npj Climate and Atmospheric Science*, 4, 1–9.

https://doi.org/10.1038/s41612-021-00189-4

Stocker, T.F. et al., 2013. Climate change 2013: The physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on
Climate Change. Cambridge: Cambridge University Press. Available from: https://www.ipcc.ch/report/ar5/wg1

Stott, P. et al., 2013. The upper end of climate model temperature projections is inconsistent with past warming. *Environmental Research Letters*, 8, 014024.
https://doi.org/10.1088/1748-9326/8/1/014024

Swanson, D.A., Tayman, J., Bryan, T.M., 2011. MAPE-R: a rescaled measure of accuracy for cross-sectional subnational population forecasts. *Journal of Population Research*, 28, 225–243. https://doi.org/10.1007/s12546-011-9054-5

Taheriyoun, A.R., Moghimbeygi, M., 2017. Visual information and expert's idea in Hurst index estimation of the fractional Brownian motion using a diffusion type approximation. *Scientific Reports*, 7, 42482. https://doi.org/10.1038/srep42482

Terzago, S., Fratianni, S., Cremonini, R., 2013. Winter precipitation in Western Italian Alps (1926–2010). *Meteorology and Atmospheric Physics*, 119, 125–136. https://doi.org/10.1007/s00703-012-0231-7

Tommasi, D. et al., 2017. Managing living marine resources in a dynamic environment: The role of seasonal to decadal climate forecasts. *Progress in Oceanography*, 152, 15–49.
https://doi.org/10.1016/j pocean.2016.12.011

Trenberth, K.E., Paolino, D.A., 1980. The Northern Hemisphere sea-level pressure data set: Trends, errors and discontinuities. *Monthly Weather Review*, 108, 855–872.
https://doi.org/10.1175/1520-0493(1980)108<0855:TNHSLP>2.0.CO;2

Valenzuela, O. et al., 2008. Hybridization of intelligent techniques and ARIMA models for time series prediction. *Fuzzy Sets and Systems, Theme: Fuzzy Models and Approximation Methods*, 159, 821–845. https://doi.org/10.1016/j.fss.2007.11.003

Valt, M., Cianfarra, P., 2010. Recent snow cover variability in the Italian Alps. *Cold Regions Science and Technology*, 64, 146-157. https://doi.org/10.1016/j.coldregions.2010.08.008
Vishwakarma, D.K. et al., 2023. Forecasting of stage-discharge in a non-perennial river using machine learning with gamma test. *Heliyon*, 9, e16290.

https://doi.org/10.1016/j.heliyon.2023.e16290

Wang, W. et al., 2005. Testing and modelling autoregressive conditional heteroskedasticity of streamflow processes. *Nonlinear Processes in Geophysics*, 12, 55–66.

https://doi.org/10.5194/npg-12-55-2005

Wang, Z.-Y., Qiu, J., Li, F.-F., 2018. Hybrid models combining EMD/EEMD and ARIMA for long-term streamflow forecasting. *Water*, 10, 853. https://doi.org/10.3390/w10070853

Xie, S.-P. et al., 2015. Towards predictive understanding of regional climate change. *Nature Climate Change*, 5, 921–930. https://doi.org/10.1038/NCLIMATE2689

Yang, Q. et al., 2023. Impact of the shrinkage of Arctic Sea ice on Eurasian snow cover changes in 1979–2021. *Advances in Atmospheric Sciences* 40, 2183–2194.

https://doi.org/10.1007/s00376-023-2272-x

Yuan, Q. et al., 2023. Assessing uncertainty in hydrological projections arising from local-scale internal variability of climate. *Journal of Hydrology*, 620, 129415.

https://doi.org/10.1016/j.jhydrol.2023.129415

Zhang, J. et al., 2023. Impacts of the Arctic stratospheric polar vortex changes on the frontogenesis over the northern middle latitudes during winter. *Atmospheric Research*, 289, 106751. https://doi.org/10.1016/j.atmosres.2023.106751
# Appendix (Teleconnection-Climate Pattern Index forecasting)

## Initial modelling rules

Extrapolating future values from past data requires a nuanced approach, distinct from the strategy employed in the NSD forecasting method discussed earlier. A comprehensive understanding of the time-series behaviour necessitates a basic visual examination of the data. This involves identifying cycles, analysing statistical patterns, and discerning trends, cycles and oscillations around a mean, as well as accounting for chaotic and unpredictable changes. To establish stationarity, a careful investigation into variations in both mean and variance, along with an assessment of heteroscedasticity in the data, becomes imperative.

Commencing an exploratory data study, we created a time-series of the Teleconnection-Climate Pattern Index (TCPI) spanning from 1884 to 2022 (Figure A1). In Figure A1a, the annual TCPI development is showcased, with the Mann-Kendall test (shown by the gray dashed linear trend) robustly suggesting a declining trend (p<0.01). The residuals of this trend are visually represented in Figure A1b.

Despite exhibiting some cyclical oscillations, the TCPI variance displays a decreasing linear trend (Figure A1c). In summary, the first and second orders of the TCPI time-series appear non-stationary. Consequently, our investigation shifts towards models capable of accommodating heteroscedasticity in the ARCH model simultaneously, ensuring autocorrelation is not constant across time. Following that, a test is employed to determine the presence of heteroscedasticity by scrutinising the residual sequence (Figure A1b).

With a McLeod-Li test yielding a p-value of 0.548, indicating the absence of heteroscedasticity, we identified the suitable PARMA(1,1)*m model for matching the TCPI data.

## Model selection
The reconstruction of the Teleconnection-Climate Pattern Index (TCPI) for the period 2022-2060 followed a distinct approach compared to that used for the NSD time-series. In this case, an autoregressive model was used without the incorporation of exogenous support information. In addition, the results of McLeod and Li (1983) test indicate that the TCPI time-series is not heteroskedastic.

The location component is reproduced by the algorithm described in Section 2.4. For the log-scale component, which does not account for heteroscedasticity, we have:

$$ \sigma(S) = \exp(\omega + \lambda + \epsilon(t)) \qquad (A1) $$

Here, the leverage term $\lambda$ accounts for the asymmetry in the data by considering how past residuals can affect future volatility beyond the baseline value $\omega$ and the residual term $\epsilon(t)$ analysed through the score term. For the log-scale component, the inclusion of the leverage effect in the EGARC model was beneficial, especially given the observed asymmetry in the data. Subsequent validation of the PARMA(TVAR)-EGARC model further confirmed its effectiveness within this simplified framework.

## Parameterisation and evaluation of the PARMA(TVAR)-EGARC model

Figure A2 depicts the observed and modeled data evolution using the PARMA[1,1]25 training period (1884-1991) and validation period (1992-2022). The closely aligned timings in the validation period closely reflect those in the training period, indicating agreement between the observation (blue) and model (orange) lines. Table A1 displays the error statistics for both the training and validation periods. The final 30 data values at the end of the time-series were utilised to compute model validation error statistics. Confirming prediction mistakes is critical for ensuring forecast accuracy. The first three statistics are used to assess the size of errors. Since the results from the validation period (1992-2022) are
comparable to the training period (1884-1991), the model is expected to perform as predicted in future projections, as shown in Table A1.

The statistics for both the training and validation stages are shown in Table A1 (box). Specifically, during the training stage, the Root Mean Square Error (RMSE) and the Mean Absolute Error (MAE) are 0.93 and 0.69, respectively, while the Mean Absolute Percentage Error (MAPE) is 8.52, with a Log Loss of 1.33.

The performance metrics (RMSE, MAE and MAPE) in the validation stage show the same level of accuracy as observed in the training step. MAPE, in particular, remains relatively low in both training and validation stages, with values around 8. A MAPE of less than 10% in the PARMAX(TVAR)-EGARC model indicates a high degree of accuracy in predictions, especially when aiming to capture long-term trends in CSD (Ostertagova and Ostertag, 2012). However, the performance metrics present values that are not entirely accurate. Specifically, Pearson's correlation coefficient (r) and KGE are below the 0.5 value, with 0.42 and 0.40, respectively.

At this stage, a more thorough analysis of the results becomes imperative. We conducted a detailed investigation using the Kolmogorov-Smirnov test to compare the distributions of the observed and predicted NSD samples during the validation stage. The calculation of the maximum distance between the cumulative distributions of the two samples yielded an estimated overall statistic (DN) of 0.27, a two-sided large sample K-S statistic of 1.03 and an approximate p-value of 0.237. The most crucial aspect is the estimated p-value, and since it is greater than 0.05, there is no statistically significant deviation between the two samples at the 95% confidence level. Consequently, we conclude that the results are suitable for practical application in the use of the PARMAX(TVAR)-EGARCH model for NSD decadal forecasting.
<figure><img src="image_12.png" /></figure>
<table><caption>**Table A1.** SARMA(TVAR)-EGARC model performance for the training (1884-1991) and validation (1992-2022) periods for the Teleconnection-Climate Pattern Index.</caption><thead><tr><th rowspan="2">Statistics</th><th colspan="2">PARMA(TVAR)-EGARC Model</th></tr><tr><th>Training period (1884-1991)</th><th>Validation period (1992-2022)</th></tr></thead><tbody><tr><td>RMSE (days)</td><td>0.93</td><td>1.17</td></tr><tr><td>MAE (days)</td><td>0.69</td><td>0.83</td></tr><tr><td>MAPE (%)</td><td>8.52</td><td>8.20</td></tr><tr><td>Log Loss</td><td>1.33</td><td>1.53</td></tr><tr><td>r</td><td>0.61</td><td>0.42</td></tr><tr><td>KGE</td><td>0.50</td><td>0.40</td></tr></tbody></table>
<table><caption>**Table 1.** Performance of the PARMAX(TVAR)-EGARCH model during the training period (1884-1991) and the validation period (1992-2022), compared with the PARMA model's performance for the Montevergine snowfall days time-series. RMSE: root mean square error (optimum, 0 to +∞, worst); MAE: mean absolute error (optimum, 0 to +∞, worst); MAPE: mean absolute percentage error (optimum, 0 to +∞, worst); Log Loss: logarithmic loss (optimum, 0 to +∞, worst); r=Pearson's correlation coefficient (optimum, +1; worst: -1); KGE: Kling-Gupta Efficiency index (worst, -∞ to 1, optimum).</caption><thead><tr><th rowspan="2">Statistics</th><th colspan="2">PARMAX(TVAR)-EGARCH Model</th><th>PARMA(TVAR) Model</th></tr><tr><th>Training period<br>(1884-1991)</th><th>Validation period<br>(1992-2022)</th><th>Validation period<br>(1992-2022)</th></tr></thead><tbody><tr><td>RMSE (days)</td><td>8.5</td><td>6.4</td><td>11.5</td></tr><tr><td>MAE (days)</td><td>6.5</td><td>4.7</td><td>9.1</td></tr><tr><td>MAPE (%)</td><td>16.1</td><td>19.1</td><td>46.4</td></tr><tr><td>Log Loss</td><td>3.4</td><td>3.3</td><td>4.2</td></tr><tr><td>r</td><td>0.77</td><td>0.68</td><td>0.16</td></tr><tr><td>KGE</td><td>0.70</td><td>0.68</td><td>0.09</td></tr></tbody></table>
## Figure captions

**Figure 1.** Schematic illustrating key features through which snowfall and snow cover influence climate, hydrological and ecological processes. Snow, a critical component of the Earth's climate system, spans approximately 46 km² of the Earth's surface, affecting regional weather patterns and playing a vital role in regulating the Earth's surface and atmosphere (adapted from Meredith *et al.* 2019, accessible for free through https://www.ipcc.ch/srocc/chapter/chapter-3-2).

**Figure 2.** Historical evolution of snow indicators. a. Winter snow severity index (SSI) reconstructed for the period 1200-2020 AD in the central Mediterranean area (based on data from Diodato *et al.*, 2019); b. Winter (October-April) snow-cover fraction evolution during the 1966-2022 period for Medieterranean area, and c. Annual mean percentage of days snow-cover during the period 1971-2010 in southern Europe. Graphs b and c were generated using data compiled from the Rutgers University Global Snow Lab, accessed through Climate Explorer (http://climate.rutgers.edu/snowcover-chart_clim.php?ui_month=8).

**Figure 3.** Environmental setting and snow outline. a. Scenic European backdrop of Italy (white square); b. Mean annual snow depth across Italy, with c. a focussed portrayal, reflecting data averaged during 1921-1960 (source: http://www.nimbus.it/liguria/rlm15/neve.htm), d. A 3-D perspective of Partenio Mountain featuring Montevergine Abbey perched near the summit, and e. Close-up of its Observatory (arranged from Gubitosa 1990). Maps in a. and d. were arranged under a free license from CESIUM (https://ion.cesium.com).
**Figure 4.** Montevergine number of snowfall days (NSD) data (1884-2022). **a.** Time-series of annual NSD; **b.** Detrended residuals; **c.** Variance estimates of NSD with an 11-year moving window.

**Figure 5.** Analysis of standardised residuals and predictability of the number of snowfall days (NSD) spanning from 1884 to 2022. **a.** Autocorrelation function (ACF); **b.** Partial autocorrelation function (PACF), featuring horizontal black lines indicating the 95% confidence level limits; **c.** Time-series attractor in the phase-space domain, showing the trajectory of NSD (blue) and distinctive ellipse-like branching (red contour).

**Figure 6.** Visual representation of the workflow illustrating the criteria for the PARMAX(TVAR)-EGARCH model criteria employed in predicting the Montevergine time-series of NSD.

**Figure 7.** Model processes in forecasting Montevergine number of snowfall days (NSD) with an exogenous variable. **a.** Scatterplot depicting the relationship between annual NSD (1884-2022) and the exogenous teleconnection climate pattern index (TCPI). The inner (deep pink area) and outer (light pink area) bounds indicate the 90% and 95% confidence limits, respectively. The periodicity scale of the signal observed on the y-axis corresponds to 2.0 per value on the x-axis; **b.** Spectral frequency density of the snowfall days time-series, highlighting the likely window (red box) where the PARMAX(TVAR)-EGARCH model cycle is located; **c.** Modelled (orange line) and observed (blue line) annual NSD time-series during the validation stage.
**Figure 8.** Residuals of the PARMAX(TVAR)-EGARCH model for the period 1884-2022. (a) Histogram and fitted normal density; (b) QQ-plot (dots) with confidence interval (blue dashed lines) of residuals.

**Figure 9.** Decadal prediction of snowfall days using the PARMAX(TVAR)-EGARCH model. Time evolution of observed (blue line, 1961-2022) and forecasted (orange line, 2023-2060) annual number of snowfall days at Montevergine Observatory, along with the corresponding long-term mean (thin gray line) and ±standard deviation (σ, blue and orange thick lines, respectively). The gray band surrounding the forecast line indicates uncertainty estimates derived from a 5-year moving window of standard deviation. Additionally, the dashed gray line represents the annual trend of the 5-year running ensemble mean of heavier snow load for Alto Adige under RCP8.5 (Faggian et al. 2024). Data on the left Y-axis are presented in log units.

**Figure A1.** Teleconnection-Climate Pattern Index (TCPI) data (1884-2022). **a.** Annual TCPI time-series with a regression line; **b.** Residuals of the TCPI regression; **c.** Variance estimates of TCPI with an 11-year moving window with regression line.

**Figure A2.** **a.** Modelled (orange line) and observed (blue line) annual Teleconnection-Climate Pattern Index (TCPI) time-series during the validation stage; **b.** Modelled (orange line) and observed (blue line) annual TCPI time-series during the forecasting stage.
<figure><img src="image_13.png" /><figcaption>Figure 1</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 2</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 3</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 6</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 7</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 8</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 9</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure A1</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure A2</figcaption></figure>