# Climate-assisted data-driven decadal snowfall predictions in the Swiss foothills

Nazzareno Diodato, Fredrik Charpentier Ljungqvist, Gianni Bellocchi

► To cite this version:

Nazzareno Diodato, Fredrik Charpentier Ljungqvist, Gianni Bellocchi. Climate-assisted data-driven decadal snowfall predictions in the Swiss foothills. PLOS Climate, 2025, 4 (7), pp.e0000592. 10.1371/journal.pclm.0000592. hal-05163270

HAL Id: hal-05163270

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

Submitted on 15 Jul 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>
RESEARCH ARTICLE

# Climate-assisted data-driven decadal snowfall predictions in the Swiss foothills

Nazzareno Diodato¹, Fredrik Charpentier Ljungqvist¹,²,³*, Gianni Bellocchi¹,⁴

1 Met European Research Observatory – International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy, 2 Department of History, Stockholm University, Stockholm, Sweden, 3 Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden, 4 Université Clermont Auvergne, INRAE, VetAgro Sup, UREP, Clermont-Ferrand, France

* fredrik.c.l@historia.su.se

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

## Abstract

Decadal-scale climate predictability is crucial for societal planning, not least in regions sensitive to winter extremes. This study predicts the number of heavy snowfall days in the Swiss Pre-Alpine Region (SPAR) through 2060, using a time-varying autoregressive model based on data from 1884 to 2023. The model integrates a pattern index that combines both large-scale (Arctic Oscillation and Dipole Model Index) and regional-scale (spring–winter temperature differential) climate forcings, capturing lagged seasonal effects. The results indicate a slight upward trend – about one to two additional heavy snowfall days by the 2050s – though not statistically significant, and consistent with regionally downscaled climate models. After 2045, variability is expected to rise, periods of snowfall deficits affecting year groups and a cluster of exceedance years emerging towards the end of the projection period. While extreme snowfall events are projected to become less frequent in the Northern Hemisphere, their intensity is unlikely to diminish. These findings enhance understanding of snowfall dynamics in the SPAR and contribute to broader insights into cryospheric changes. Data-driven models such as this one are valuable tools to contextualise historical hydroclimatic drivers and assessing inter-annual and inter-decadal variability in regional climate projections and their societal implications.

## 1. Introduction

Snow has important ecological, social and economic impacts across various regions of the world [1,2]. Heavy snowfall events are expected to have growing impact in the next decades [3–5]. These events interact with the global water cycle at different time-scales [6,7], influencing surface temperatures and the hydrological cycle, particularly in snow-dominated mountain regions [8, Fig 1]. Consequently, understanding and predicting snow variability is essential for accurate climate projections, especially in regions like the European Alps that are particularly vulnerable to climate

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

**Citation:** Diodato N, Ljungqvist FC, Bellocchi G (2025) Climate-assisted data-driven decadal snowfall predictions in the Swiss foothills. PLOS Clim 4(7): e0000592. https://doi.org/10.1371/journal.pclm.0000592

**Editor:** Sher Muhammad, ICIMOD: International Centre for Integrated Mountain Development, NEPAL

**Received:** January 30, 2025

**Accepted:** June 2, 2025

**Published:** July 2, 2025

**Copyright:** © 2025 Diodato et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** All data used in this study are freely available. The graphs produced are original creations of the authors from the original sources. The full set of data supporting the conclusions of this study (observed and forecasted snowfall data), is provided in the supplementary file to this article (S1 Table).
<figure><img src="image_4.png" /></figure>

**Funding:** This research was performed as an investigator-driven study without financial support for N. D. and G. B. F. C. L. was supported by the Swedish Research Council (Vetenskapsrådet, grant no. 2023-00605), the Marianne and Marcus Wallenberg Foundation (grant no. MMW 2022-0114), and the Centre for Advanced Study (CAS) at the Norwegian Academy of Science and Letters, which funded and hosted the research project "The Nordic Little Ice Age" during the 2024/25 academic year. Open access publication funding for this article was provided by Stockholm University.

**Competing interests:** The authors have declared that no competing interests exist.

variability due to their sensitivity to temperature shifts, dependence on snowpack for water resources, and exposure to associated ecological and socio-economic impacts [9–12]. However, despite the capabilities of advanced climate models such as Global Circulation Models (GCMs) and Regional Climate Models (RCMs), uncertainties persist. While RCMs offer improved spatial resolution by incorporating regional topography [13] and may be valuable in capturing regional climate features, comparing their performance with GCMs in simulating specific phenomena, such as snowfall events, remains challenging.

Studies like Fereday et al. [14] highlight the influence of atmospheric dynamics as a major source of uncertainty in climate predictions. While physically based-models account for complex interactions between variables, these interactions can introduce uncertainties, particularly when the underlying processes are not fully understood or when input data are limited [15]. These challenges highlight the need for alternative approaches, such as statistical models, to address snowfall prediction. Statistical models, particularly data-driven models (DDMs), offer a complementary way to explore climate variability and predict future trends. Focusing on teleconnection patterns, which describe large-scale atmospheric phenomena influencing regional climate, DDMs can capture internal climate variability and provide insights into snowfall predictions [16,17].

While existing literature has explored snowfall patterns in various regions, including the European Alps [21–23], and the relationships between snowfall and climate patterns [24,25], there is a gap in combining teleconnection patterns with multiple climate variables to predict specific aspects like snowfall days [26,27] or snow cover duration [28]. This gap exists partly because GCMs and RCMs often struggle with accurately capturing snowfall variability due to computational limitations [29,30].

Understanding climate patterns has long been essential for guiding land use decisions, particularly in agriculture. The principle is not new: classical writers like the Roman poet Virgil (70–19 BCE) emphasised the value of anticipating climatic variability. In *Georgics*, Virgil advised farmers to observe wind patterns, past harvests and the character of the land before sowing crops (Solari, [31], pp. 79–80, translation from https://www.poetryintranslation.com/PITBR/Latin/VirgilGeorgicsI.php):

> [...] know the winds, and the varying mood of the sky, and note our native fields,
> and the qualities of the place, and what each region grows and what it rejects [...].

This ancient perspective resonates with the modern aim of this study – to model regional climatic behaviour using historical data and climatic indices in order to guide future preparedness. Just as Virgil highlighted the role of environmental knowledge in planning, our study reinforces the need for climate-informed modelling to support decision-making in snow-sensitive regions such as the European Alps. DDMs, which rely on past data to simulate natural climate variability, offer a promising alternative for predicting snow trends over decadal time-frames [32]. As Chatfield [33] noted, model selection is both an art and a science. Mulligan and Wainwright [34] similarly described modelling as an art that requires expertise, intuition and mathematical
<figure><img src="image_5.png" /><figcaption>**Fig 1.** Schematic picture depicting some of the fundamental aspects and processes by which the cryosphere and atmosphere in snow-covered regions impact global climate, ecological, and social systems (modified from Intergovernmental Panel on Climate Change [18]; https://www.ipcc.ch/srocc/chapter/chapter-3-2). Variations in snowfall and snow-cover, mediated by complex teleconnections and feedbacks, have significant implications for various Earth system components, including water supplies, ecosystems and global temperature dynamics [19]. Teleconnected weather patterns provide critical insights into the drivers behind changes in snowfall, both encompassing long-term trends and short-term atmospheric circulation variability [20], thereby enhancing climate prediction capabilities.</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g001

creativity. This is particularly evident in DDMs, where identifying meaningful interactions between variables requires not only analytical skill but also interpretive insight.

Snowfall, as a significant component of the Earth's energy budget, exerts a distinct forcing on the climate system through the albedo effect, latent and sensible heat fluxes and its role in modulating atmospheric circulation patterns [35]. Variations in snowfall from year to year, driven by factors such as temperature, precipitation and large-scale atmospheric circulation [36], can influence regional and global climate trends through these mechanisms [37]. While recent research, primarily leveraging quantitative methods and large-scale datasets, has linked past weather patterns to a wide array of social and economic outcomes, providing compelling evidence of the significant impacts of climate on human systems [38], there remains a notable gap in research focusing on the long-term variability of snowfall, particularly heavy snowfall events.

By providing extensive, uninterrupted estimates of unforced natural variability, historical records are invaluable for understanding natural climate variability and the influence of natural forcing on climate dynamics. DDMs, which involve the training of autoregressive models based on observed time-series data, can be employed to capture this internal variability, allowing for some extrapolation into future climate scenarios [39]. Recent advancements in statistical modelling, such as Autoregressive Integrated Moving Average (ARIMA) models, have been applied to precipitation forecasting [40–42], but snowfall predictions remain relatively understudied [43]. Snow climate models often demonstrate proficiency in forecasting winter temperature and precipitation patterns. However, they frequently exhibit limitations in accurately forecasting the year-to-year variability of snowfall, including both total snowfall amounts and spatial variations in the magnitude and
<figure><img src="image_6.png" /></figure>

extent of snow cover [44]. Studies like Clark et al. [45] and McCray et al. [46] have explored specific aspects of snow-
fall variability in areas of eastern North America. Hammer et al. [47] provided seasonal snow depth forecasts in Norway
applying a Seasonal Autoregressive Moving Average (SARMA) model, and another notable advancement from Koehler
et al. [12] combines different statistical approaches, including autoregressive models, for decadal snow-cover line fore-
casting in the European Alps between 2022 and 2029. These mountains have experienced a marked decrease in snowfall
days, snow depth, and snow cover since the 1960s, particularly at low and intermediate elevations [44]. This decline has
been observed at many meteorological stations in the Swiss region compared to the reference period from 1960 to 1985
[48]. Consequently, shorter winter periods and a shift from snow-dominated to rain-dominated conditions have ensued
[21], impacting regions differently. This transformation has led to increased winter discharge [49,50] and erosion [51]. The
regional and seasonal nuances of these changes underscore the complexity and varied impacts of shifting snowfall pat-
terns in different parts of the world.

This study addresses the gap in snowfall prediction by focusing on the number of heavy snowfall days (NHSD), defined as the number of days between December and March with a snow depth equal or greater than 0.05 m. This metric is based on data for the Swiss Pre-Alpine Region (SPAR), originally derived by Marty [21] and updated by Diodato et al. [52]. Using a time-varying (TVAR) periodic autoregressive moving average (PARMA) approach with exogenous (X) climate forcing and conditional standard deviation (CSD), the study aims to project NHSD trends until 2060. By incorporating both large-scale (Arctic Oscillation and Dipole Model Index) and small-scale (temperature shifts) factors, the model has a better chance of capturing the complex interplay of climate variables. The use of time-varying frequency analysis enhances the ability of the model to accurately forecast the NHSD under chaotic conditions [53]. In addition, exogenous control can mitigate non-linearity and eliminate shifts in the autoregressive model, leading to more reliable and accurate predictions [54]. To develop a PARMAX(TVAR)-CSD model, the study uses a reliable long-term time-series (1836–2023) of NHSD for SPAR [21,52]. Due to the availability of relevant climate forcing data, the analysis focuses on the period from 1884 onwards.

## 2. Data and methods

### 2.1. Environmental setting

The Swiss Pre-Alpine Region (SPAR) is located on the northern slope of the Alps in central-western Europe (Fig 2a, b).
This subalpine strip, with elevations ranging from approximately 800–1500 m a.s.l., experiences annual snowfall lasting
between 20 and 80 days (Fig 2c), resulting in a consistent winter snow period from November to April (Fig 2d).

SPAR's landscape has been shaped by centuries of human activities, such as alpine grazing and deforestation, with
archaeological evidence tracing back to the Bronze Age (ca. 3300–1200 BCE). However, compared to regions with inten-
sive farming or dense settlements, SPAR has retained a relatively limited extent of landscape alteration [55]. Three distinct
physical regions characterise the area: (1) alpine region to the south, constituting over 60% of the total area; (2) pre-alpine
central region below 800 m a.s.l., occupying 27% of the area (Fig 2b, yellow zones); (3) northern region, more extensively
influenced by agriculture and urbanisation.

The pre-alpine central region, nestled within the Swiss Plateau, displays a continuous landscape of hills, steep valleys and small plains. Although the SPAR covers a relatively small area, it exhibits substantial climatic variability driven by local features such as large lakes and diverse slope exposures, both moderating snowfall patterns [56]. The SPAR typically encounters heavy snowfall (more than two days per week) from December to February, which accounts for a substantial proportion (likely in the range of 60% to 80%) of the total snow amount and helps maintain a sustained snow cover into May (Fig 2d, black line and white bars).

Snowfall patterns follow a northwest-southeast gradient, with maxima occurring in areas exposed to cold and moist north-westerly flows [22], with similarities between the French and Swiss Alps. In addition, local factors such as geographical characteristics and atmospheric convection contribute to snowfall variability [57]. For example, blocking high-pressure systems near Greenland can channel cold air into Western Europe, leading to increased snowfall in the Alps [58].
<figure><img src="image_7.png" /><figcaption>**Fig 2. Environmental setting and snow outline.** a. Swiss region setting (red square) adapted from Shaded Relief (http://www.shadedrelief.com/natural2/globes/africa.jpg); b. Pre-alpine region (yellow areas) adapted from the National Centre for Climate Services, MeteoSwiss (https://www.nccs.admin.ch/nccs/en/home/regions/grossregionen/pre-alps.html); c. Mean annual snowfall days across Switzerland (1981–2010), with the Pre-Alpine region outlined in white, adapted from the National Centre for Climate Services-MeteoSwiss (https://www.nccs.admin.ch/nccs/en/home/climate-change-and-impacts/swiss-climate-change-scenarios/facts-and-figures/climate-iindicators.html); d. Winter cycle of heavy snowfall days (black line), snow amount (histogram) and rainfall days (dashed blue line) per weeks in each month from October to May, averaged on the period 2007–2022 for the Swiss Valley, adapted from Snow-Forecast (https://www.snow-forecast.com/resorts/Swiss-Valley/history).</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g002

## 2.2. Explanatory data analysis and predictability

Given SPAR's complex environmental context, understanding regional snowfall patterns is crucial for assessing temporal changes. Thus, we applied exploratory data analysis techniques to examine snow variability, an approach that emphasises graphical data visualisation over purely statistical analyses [59]. Exploratory data analysis facilitates pattern and anomaly detection in datasets through visual tools like graphs and plots. The climate system exhibits variability on different time-scales, including precipitation patterns. However, the structure of this variability is not random and can be characterised by scaling relationships, such as power laws observed in probability density distributions and autocorrelation functions [60].

Our analysis investigated NHSD at various time-scales, capturing their periodicity and autocorrelation features. Initial exploration of the data revealed annual fluctuations, with the Mann-Kendall test (p<0.01) and the Kendall-tau test
<figure><img src="image_8.png" /></figure>

(p < 0.01) indicating a significant declining trend in NHSD (Fig 3a). The Mann-Kendall and Kendall-tau tests are non-parametric methods used to detect trends in time-series data without assuming a specific distribution. Further analysis showed interannual variability, where residuals exhibited oscillatory behaviour captured through an 88-year cosinusoidal filter (Fig 3b). A cosinusoidal filter is a smoothing technique used to extract long-term cyclical components from time-series data. It reduces short-term fluctuations and highlights oscillatory patterns, making it particularly well-suited for identifying interannual variability in the NHSD dataset [61]. This oscillatory behaviour was also evident in the conditional variance estimates, which displayed periodic components across both short-term (~40 years) and long-term time-scales (Fig 3c). Conditional variance estimates how variability evolves in a time-series. Using an 11-year moving window, this measure helps uncover temporal shifts in fluctuation intensity, particularly in non-stationary data like NHSD. These periodicities suggest suitability for predictive modelling using the PARMAX framework for these time-scales.

To model the time-series dynamics of NHSD, we used a suite of advanced techniques – including the autocorrelation function (ACF), partial autocorrelation function (PACF), Hurst exponent analysis and attractor dynamics – which collectively provide robust support for understanding autocorrelation and long-term dependencies [62,63]. The autocorrelation function shows how current values in a time-series are related to past values, helping to detect persistence or repetition. The partial autocorrelation function quantifies the direct influence of prior values on the current value in a time-series, excluding the effects of intervening lags. It is particularly useful for identifying the true order of autoregressive processes. The Hurst exponent measures the long-term memory of a time-series, indicating whether trends are likely to persist or reverse over time. Attractor dynamics describe how a system evolves over time by mapping its behaviour in phase-space, a geometric representation that reveals stability, cycles and recurrence. In this space, patterns such as ellipses suggest deterministic and potentially predictable system dynamics. Together, these methods offer comprehensive insights into the periodicity, persistence and predictive potential of the NHSD time-series. The ACF and PACF are indispensable tools for identifying the nature of dependencies in time-series data. The ACF quantifies the correlation between the NHSD time-series and its lagged versions, enabling the detection of persistent patterns (Fig 4a).

<figure><img src="image_9.png" /><figcaption>**Fig 3.** Swiss pre-Alpine region: heavy snowfall days (NHSD) data (1884–2023). a. Annual NHSD (blue line) with a significant linear trend (grey dashed line); b. Detrended residuals (violet line) overlaid with an 88-year cosinusoidal filter (red line), emphasising periodic trends; c. NHSD conditional variance, estimated using an 11-year moving window. The vertical scale of the band-pass filter in b is magnified threefold to improve phase visibility in relation to periodic trends.</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g003
<figure><img src="image_10.png" /><figcaption>**Fig 4.** Analysis of the number of heavy snowfall days across the Swiss pre-Alpine region. a. Autocorrelation function (ACF), and b. Partial autocorrelation function (PACF) on the standardised residuals; featuring horizontal black lines indicating the 95% confidence level limits; c. Time-series attractor in the phase-space domain, showing the trajectory of NHSD (blue) and the distinctive ellipse-like branching (red contour) for the period 1884–2023.</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g004

Notably, a significant spike in the ACF at lag $k=1$ highlights a correlation between consecutive years' NHSD values. The alignment of the short-memory behaviour observed in the autocorrelation with the long-term persistence measured by the Hurst exponent ($H>0.5$) underscores the reliability of the NHSD's temporal dynamics [64]. The ACF displays notable correlations at several lags, indicative of underlying dependencies in the NHSD time-series, suggesting that past NHSD values influence future dynamics. The PACF complements the ACF by isolating direct relationships between NHSD values at specific lags, controlling for intermediate lags' influence (Fig 4b). This distinction is critical, as indirect correlations in the ACF may obscure the true dynamics. Analysis of the PACF revealed irregular patterns, pointing to the suitability of an ARMA($P,Q$) structure for capturing the time-series dependencies comprehensively [65]. The identified irregular patterns in the ACF and PACF informed the selection of a PARMA(1,1) model with a monthly cycle ($m=41$). This choice aligns with the need to accommodate periodic and exogenous influences, such as climate predictors, that drive the NHSD dynamics. The incorporation of these external factors reflects a sophisticated understanding of the interactions between intrinsic time-series behaviour and external climate variability. The model's robustness was further confirmed by the McLeod and Li [66], which indicates the absence of heteroscedasticity ($p=0.80$). The McLeod and Li test checks whether the variability in a time-series is consistent over time (homoscedasticity).

The Hurst exponent $H$ was calculated using the Rescaled Range Analysis (R/S) method [67], yielding $H=0.77$. Rescaled Range Analysis is used to estimate the Hurst exponent and assess whether a time-series shows random or persistent behaviour. This value signifies strong long-term dependency, reflecting a time-series with persistent behaviour where historical trends are likely to continue [68]. The relationship between $H$ and the ACF is particularly noteworthy: persistent time-series often display an ACF with slow decay, a characteristic of long memory processes [69]. This observation supports the use of fractal-based methods in analysing NHSD dynamics and further justifies the choice of models designed to accommodate such persistence, particularly those that incorporate periodicity and exogenous climate drivers.

Phase-space attractor analysis revealed an ellipse-like trajectory for NHSD (Fig 4c), a pattern indicative of structured and potentially predictable dynamics. The presence of such an attractor suggests resilience to disturbances and emphasises the system's inherent regularities. This structured behaviour reinforces the presence of deterministic components within the NHSD time-series, aligning with the fractal characteristics suggested by the Hurst exponent and ACF/PACF
analysis. The attractor's stability underlines the suitability of the PARMAX(TVAR)-CSD framework for forecasting NHSD patterns. By capturing both short-term oscillations and long-term trends, this model can effectively simulate the NHSD time-series behaviour.

## 2.3. PARMAX(TVAR)-CSD model framework

The PARMAX(TVAR)-CSD model was selected for NHSD trend prediction due to its robustness in modelling time-series data with exogenous climatic influences [34]. Fig 5 outlines the model workflow, highlighting key processing steps.

The random variable $y_t$, representing NHSD at time $t$, is modelled using both its past values and a deterministic seasonal component with fixed period $m$. To incorporate exogenous variables, the model adopts an additive structure, consisting of two main components: a stochastic autoregressive component and a deterministic seasonal cycle.

The location component is given by [70]:

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

Here, $y_t$ is the NHSD value at time $t$ (the step size for $t$ is one year.); $S_{t-i}(ARX)$ denotes the random autoregressive (AR) process with exogenous inputs (ARX), capturing short-term dependencies and external climatic effects; $m_{t-jm}$ is the deterministic, periodic component with period $m$, representing seasonal trends ($m=41$ years corresponds to the periodicity observed in the dataset).

<figure><img src="image_11.png" /><figcaption>Fig 5. Visual representation of the workflow illustrating the criteria for the PARMAX(TVAR)-CSD model criteria employed in predicting the NHSD for the Swiss pre-Alpine region (image of snowy landscape derived from Freepik, https://www.freepik.com/free-vector/watercolor-background-with-winter-landscape_1433112.htm#query=snowfall%20landscapes&position=1&from_view=search&track=ais&uuid=3a748666-b7de-4ee9-b1ba-1bd3383c2928).</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g005
To further refine the model, we introduced an expanded expression for the location component:

$$ \mu(L) = \omega + AR2(1) \cdot \Phi_{1}^{(t-1)} + AR2(2) \cdot \Phi_{2}^{(t-2)} + P(m) \cdot \kappa + X(t) \cdot \beta1 + \varepsilon(t) \quad (2) $$

In this formulation, μ(L) is the expected value (mean) of the NHSD time-series at lag L; ω is a constant term representing the baseline value of the time-series when all other factors are zero; AR₂(1) and AR₂(2) are the first and second lagged values (e.g., yₜ-₁ and yₜ-₂) in a second-order autoregressive process; Φ₁ = 1.0362 and Φ₂ = -0.1488 are the corresponding autoregressive coefficients, which determine how past values influence the current prediction; P(m) is a periodic function with a fixed period m = 41 years, capturing the long-term oscillatory behaviour in the data; κ = 2.000·10⁻⁸ is a scaling coefficient for the periodic component; X(t) is the exogenous input at time t, with coefficient β₁ = -11.069; ε(t) is the residual error term, representing the unexplained variance between the actual and predicted values after accounting for autoregressive, seasonal and exogenous effects. The parameter values Φ₁, Φ₂, κ and β₁ were obtained through model training using historical NHSD and climate data specific to the SPAR dataset.

The score-driven methodology offers three major advantages [71]: (i) the model's likelihood is expressed in closed form, facilitating estimation; (ii) the filtered estimates of time-varying parameters are optimal in terms of Kullback-Leibler; (iii) the model's predictive capacity is similar to parameter-driven models [72].

For the log-scale component, we used the following equation:

$$ \sigma(S) = \exp(X(t) \cdot \beta2 + \varepsilon(t)) \quad (3) $$

where σ is the standard deviation; X is the exogenous variable. β₁ in Eq. (2) and β₂ in Eq. (3) are the associated parameters. Further details and numerical solutions can be found in Basawa et al. [73], who introduced autoregressive models with periodically varying parameters to enhance the representation of long-term seasonal structures. In addition, and Lit et al. [74] outlined the implementation of score-driven models, which facilitate efficient estimation and forecasting through dynamic parameter updates derived from likelihood scores.

## 2.4. Exogenous climatic input model design

Decadal forecasts benefit from incorporating pronounced long-term climatic oscillation signals into AR models [75]. Studies have shown that predictions based on teleconnection patterns outperform those lacking such forcing [76]. This highlights the combined influence of internal climatic variability and external forcing on decadal forecasts, which is the framework used in this study. While large-scale atmospheric circulation drives heavy snowfalls, local factors like topography and thermodynamics can significantly amplify these effects [77,78]. Examples include the European Alps, the Great Lake region of North America and the islands of Japan, where local characteristics enhance precipitation and contribute to high snowfall events [79].

To improve the predictability of NHSD, we incorporated a Climate Forcing Index (CFI) into the PARMAX(TVAR)-CSD model. The CFI synthesises three critical atmospheric influences, spanning both large-scale and regional dynamics: the Arctic Oscillation (AO), the Dipole Mode Index (DMI) and the spring-winter temperature differential (TDsw). The index is defined as:

$$ CFI = AO + DMI - 0.3 \cdot TDsw + 5 \quad (4) $$

These three variables were identified following extensive testing of other candidates, such as El Niño Southern Oscillation (ENSO), the North Atlantic Oscillation (NAO), and local humidity and wind metrics. None of these alternatives consistently enhanced predictability. The empirically determined CFI, with coefficients of 0.3 and 5, showed a stronger correlation with NHSD than individual exogenous variables (r = 0.67; p < 0.05). This highlights the integrated impact of tropical sea surface
<figure><img src="image_12.png" /></figure>

temperatures (SSTs), high-latitude circulation patterns and regional temperature dynamics on snowfall extremes, aligning
with the influence of large-scale atmospheric systems over the Atlantic-Eurasia region [80].

The AO, a well-established driver of Europe's winter climate, including Switzerland's winter climate, affects the strength and positioning of westerly winds, thereby influencing moisture transport and snowfall patterns. The DMI quantifies the east-west SST gradient across the equatorial Indian Ocean, serving as a key metric of the Indian Ocean Dipole (IOD). Although primarily associated with tropical weather systems, the DMI has been linked to teleconnections that influence Eurasian and Alpine winter climates, through modulation of mid-latitude wave propagation [81]. The TDsw, a regional variable, captures localised temperature gradient variability, which strongly affects snowfall thresholds via internal atmospheric dynamics [82]. The negative weight assigned to TDsw captures its inverse relationship with snowfall, while the constant value of 5 ensures the index remains within a positive range for interpretability.

Leveraging these insights, we used the PARMAX(TVAR)-CSD to assess the combined effects of large-scale forcing (AO and DMI) and regional-scale forcing (temperature jump) through the CFI. The index was incorporated into both the location and log-scale components of the model – Eqs. (2) and (3), respectively – as the exogenous input X(t). Future CFI values were projected using a PARMA(2,1) model, where the 25-year periodicity (m=25) aligns with observed oscillatory modes in the AO. The model structure was selected via heteroskedasticity testing (p=0.911) and showed no significant distributional divergence from observations during validation (Kolmogorov-Smirnov p=0.80). Performance metrics (e.g., root mean square error of 0.56 and mean absolute error of 0.45) confirmed the robustness of the forecasting approach. The CFI was prepared in a sub-routine and entered into the information flow of NHSD forecasting (full technical details in Supporting Information).

**2.5. Model assessment**

In any prediction task, particularly one involving complex real-world variables, ensuring reliable forecasts and effective model evaluation is crucial [83]. This is because the inherent complexity often makes it challenging to consistently predict future values. Validating forecasting is thus imperative to allow users to grasp potential discrepancies between forecasts and actual observations.

For this study, the entire dataset spanning 1884–2023 was divided into three periods: training (1884–1993), validation (1994–2023) and forecast (2024–2060). The training period encompasses over three-quarters of the total data (110 out of 140 years) to ensure sufficient data for training the selected model type. This resulted in a validation period of 30 years. While a longer validation period would be ideal, limitations in historical data availability necessitated this allocation. To avoid under-representing potential anthropogenic forcings in future projections, our forecast was performed using the entire training period of the time-series, spanning 1884–2023. While the earlier years (1884–1993) experienced a lower influence from anthropogenic climate change, including the full period helps account for more recent climatic shifts and better captures the evolving influence of human activities on snowfall dynamics.

Predictions were made using the Time Series Lab [74], Score Edition software V. 1.5 (https://timeserieslab.com). Addi-
tional software tools were used for various data processing: Visual Recurrence Analysis (https://visual-recurrence-anal-
ysis.software.informer.com/4.9) and STATGRAPHICS (http://www.statpoint.net) for a variety of statistical analyses, and
CurveExpert Professional 1.6 (https://www.curveexpert.net) for creating graphical visualisation. Finally, SELFIS self-
similarity analysis software (http://www.cs.ucr.edu/tkarag/Selfis/Selfis.html) was used to calculate the Hurst exponent.

**3. Results and discussion**

**3.1. The implications of exogenous (climate) variables**

The linear relationship between the CFI and yearly NHSD is displayed in Fig 6a. Notably, most data points fall within the 95% prediction bounds, while only three data points fall outside these limits. The *p*-value in the analysis of variance (ANOVA) is less than 0.05, with Pearson's correlation coefficient equal to 0.67, suggesting a statistically significant correlation between the NHSD and the climate forcing index. Analysing each single factor of the climate forcing index at a
<figure><img src="image_13.png" /><figcaption>**Fig 6. Climate drivers of heavy snowfall in Switzerland:** a. Scatter-plot illustrating the linear relationship between the Climate Forcing Index (CFI) and the annual number of heavy snowfall days (NHSD) in the pre-Alpine region (1884–2023), with 90% (dark pink) and 95% (light pink) confidence limits; b. Spatial correlation between Arctic Oscillation (AO) sea-level pressure (https://www.atmos.colostate.edu/~davet/ao/Data/ao_index.html) and snow cover (1966–2022); c. Same as b, but for the Dipole Model Index (DMI) HadSST1 sea-surface temperature dataset (Hadley Centre Sea Surface Temperature dataset version 1; https://www.metoffice.gov.uk/hadobs/hadisst). Panels b and c are based on Rutgers University Climate Lab snow-cover data (https://climate.rutgers.edu/snowcover) and were generated using Climate Explorer (http://climexp.knmi.nl; [84]). The black box in both maps highlights the Swiss region.</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g006

time, we obtain that AO and DMI are anti-correlated with NHSD, with a Pearson's correlation coefficient equal to -0.43 and -0.32, respectively. Since the *p*-value in the tests residuals of Durbin-Watson (DW) statistic is greater than 0.05, there is no indication of serial autocorrelation in the residuals at the 95% confidence level. In determining whether the whole the terms of CFI equation can be accounted significant, notice that the highest *p*-value on the independent variables is 0.056 belonging to DMI. The *p*-value is just above 0.05, indicating that the term may not be statistically significant. However, given that the climate experience more changes and oscillations during the phase of reconstruction compared to the period used in this test, we assume a higher correlation between DMI and NHSD over a longer historical timeline. This indicates that we should not remove any variables from the CFI equation.

The spatial correlations between AO, DMI and snow cover (Fig 6b and 6c) suggest that AO and DMI influence snowfall variability across central Europe through teleconnections. A similar correlation was confirmed at the European scale for the period 1972–2006, showing that snow-cover duration is associated with large scale atmospheric patterns, particularly AO during winter [85]. The impact of DMI (IOD) was also evaluated due to its significant influence on Eurasian snow-cover duration [86]. TDsw is, instead, positively correlated, with a Pearson coefficient of *r*=0.53. As a local variable, TDsw shows the highest correlation, contrasting with AO and DMI, which represent large-scale teleconnection drivers. This correlation coefficient is comparable to the findings of Scherrer and Appenzeller [87], who identified a relationship between variations in the first snow day pattern and temperature anomalies in the Swiss Alpine region.

To assess the significance of the model terms (AO, DMI and TDsw), we merged the training and validation datasets, covering a longer time-series that spans climate conditions both before and after the recent (post-1990) global warming. The highest *p*-value among the independent variables (*p*=0.0499), corresponding to DMI, remains statistical significance (*p*<0.05). Consequently, it is advisable to retain all variables within the exogenous information model. This suggests that each term holds statistical relevance, supporting the operational use of CFI in forecasting snowfall days.

## 3.2. Accuracy and forecast assessment

The temporal behaviour during the validation period closely mirrors that of the training period, as shown by the strong agreement between the observed (blue) and modelled (orange) lines in Fig 7a. A closer inspection reveals a subtle yet
<figure><img src="image_14.png" /><figcaption>Fig 7. Comparison of modelled (orange) – Eq. (2) – and observed (blue) NHSD annual time-series in the Swiss Pre-Alpine region, with change-points in 1988 (observed) and 1986 (joint modelled).</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g007

potentially important shift: the model's local maxima tend to lag one to two years behind the observed peaks, particularly toward the end of the training period and into validation period. This misalignment may reflect a delayed model response to climate drivers not fully captured by the input variables, possibly due to lagged effects associated with atmospheric circulation patterns or seasonal temporal shifts. However, the Mann-Whitney-Pettitt homogeneity test [88] further supports the overall alignment, identifying change-points in 1988 for the observed data and 1986 for the joint modelled time-series – coinciding with the well-documented warming across Europe in the late 1980s (e.g., [89,90]).

Model performance was further validated using the last 30 data points, with error metrics (Table 1) demonstrating strong performance of the PARMAX(TVAR)-CSD model across RMSE, MAE, MAPE, r and KGE. The small discrepancies between observed and estimated values suggest limited uncertainty [91], likely due to narrow parameter sensitivity [92], as previously discussed [93,94]. In light of the low error rates, extensive sensitivity analysis may be unnecessary [95], consistent with the guidance of Sun et al. [96], and Hyndman and Athanasopoulos [97].

Although the PARMAX(TVAR)-CSD model achieved a MAPE of less than 25% (acceptable according to Swanson et al. [98]) during training, its performance declined during the validation phase, suggesting the need for further analysis. We thus performed the Kolmogorov-Smirnov (K-S) test, which showed no significant difference between observed and predicted NHSD distributions (p=0.13), reinforcing the model's reliability.

Comparing the PARMAX(TVAR)-CSD model to the simpler PARMA(TVAR) model, which lacks exogenous variables, reveals that the latter had approximately twice the RMSE and MAPE values, and showed no correlation with observations (r~0.00), highlighting the advantage of including exogenous inputs in NHSD forecasting.

# 4. Outlook and perspectives

## 4.1. Decadal forecasting of snowfall dynamics

Understanding and predicting snowfall dynamics involves considering a wide range of time-scales due to the complexity of meteorological and climatic systems. These scales span short-term fluctuations to long-term trends, making the comprehension of the NHSD framework a major challenge. However, projections show that future NHSD extremes will likely remain low, staying within the mean snowfall deficit range until around 2045, followed by remarkable interannual and
<table><caption>Table 1. Performance comparison of PARMAX(TVAR)-CSD and PARMA(TVAR)-CSD models during training (1884–1993) and validation (1994–2023) periods in the Swiss Pre-Alpine region. Sim: simulated values; Obs: observed values; RMSE: root mean square error (optimal, 0; worst: +∞); MAE: mean absolute error (optimal, 0; worst: +∞); MAPE: mean absolute percentage error (optimal, 0; worst: +∞); r = Pearson's correlation coefficient (optimum, +1; worst: -1); KGE: Kling-Gupta Efficiency index (worst, -∞ to 1, optimum); Log Loss: logarithmic loss (optimum, 0; worst: +∞).</caption><thead><tr><th rowspan="3">Statistics</th><th colspan="4">PARMAX(TVAR)-CSD Model</th><th>PARMA(TVAR)-CSD Model</th></tr><tr><th colspan="2">Training period</th><th colspan="2">Validation period</th><th rowspan="2">Validation period<br>(1994–2023)</th></tr><tr><th colspan="2">(1884–1993)</th><th colspan="2">(1994–2023)</th></tr></thead><tbody><tr><td>Mean</td><td>Sim</td><td>Obs</td><td>Sim</td><td>Obs</td><td>Sim</td></tr><tr><td></td><td>34</td><td>35</td><td>26</td><td>22</td><td>31</td></tr><tr><td>RMSE (d yr<sup>-1</sup>)</td><td colspan="2">8.1</td><td colspan="2">12.2</td><td>19.3</td></tr><tr><td>MAE (d yr<sup>-1</sup>)</td><td colspan="2">6.0</td><td colspan="2">10.0</td><td>13.7</td></tr><tr><td>MAPE (%)</td><td colspan="2">21</td><td colspan="2">110</td><td>210</td></tr><tr><td>Log Loss</td><td colspan="2">3.37</td><td colspan="2">4.17</td><td>5.84</td></tr><tr><td>r</td><td colspan="2">0.78</td><td colspan="2">0.55</td><td>~0.00</td></tr><tr><td>KGE</td><td colspan="2">0.78</td><td colspan="2">0.44</td><td>-0.13</td></tr></tbody></table>

https://doi.org/10.1371/journal.pclm.0000592.t001

intradecadal fluctuations (Fig 8, red line). Despite reaching its lowest levels between 2015 and 2035, the number of snowy days shows a slight increase, although this rise is not statistically significant (Mann-Kendall test, *p*=0.90; Kendall tau test, *p*=0.40). By the end of the forecast period, some years are projected to exceed the mean and approach the critical

<figure><img src="image_15.png" /><figcaption>**Fig 8.** Decadal prediction of heavy snowfall days. Observed data (1950–2023, blue) and forecasted values (2024–2060, red) are shown for the Swiss Pre-Alpine Region alongside the long-term mean (grey line), with the 10th and 90th percentiles (indicating snowfall deficits and heavy snowfall exceedances, respectively) based on the 1970–2023 reference period. The grey shaded band represents uncertainty, calculated using a five-year moving standard deviation. The dashed grey line indicates a five-year running ensemble mean of heavier snow loads for the northeast Italian Alps under the RCP8.5 scenario (12 km resolution; [99]).</figcaption></figure>

https://doi.org/10.1371/journal.pclm.0000592.g008
<figure><img src="image_16.png" /></figure>

threshold of heavy snowfall days – defined as the 90th percentile of NHSD observed between 1970 and 2023 – though
these levels still fall short of the historical peaks recorded prior to the 1990s.

The forecast includes 0.95 confidence intervals estimated through bootstrapping, providing robust estimates that account for non-normality. A grey band around the forecasted red line (Fig 8) represents uncertainty derived from a 5-year moving window of standard deviation, illustrating limited variability in the projections. For SPAR, an increase of NHSD is anticipated in both the higher alpine region and the northeastern Italian Alps on interdecadal and longer-term scales (Fig 8, dashed grey line). These projections align with Euro-CORDEX models [99] under the RCP8.5 emission pathway. While updated snowfall projections exist for other RCPs, incorporating them into future research can provide a more nuanced understanding of various emission scenario's impacts on snowfall patterns.

Although RCP8.5 (characterised by extremely high fossil fuel emissions) represents a less likely [100] and increas-
ingly implausible [101] emissions pathway, it serves as a valuable counter-scenario for assessing extreme situations
(e.g., [102]). Based on simulations, the northeastern Alpine region is expected to experience increased wet snow loads
(WSL) until around 2060, highlighting the need for adequate infrastructure design in high elevations, as suggested by
EURO-CORDEX regional climate model (GCM-RCM) experiments for the French Alps [103].

The potential link between major snowfall recovery in the Swiss Alps and increased mid-latitude cold extremes during recent Arctic amplification [104] remains a subject of ongoing investigation (e.g., [105]). While there is some evidence of rising cold extremes in certain regions, like the eastern United States during the 2010s, long-term data do not confirm a widespread increase at mid-latitudes. Further research is needed to fully understand these relationships. When assessing Arctic amplification's impacts on Northern Hemisphere snow, increased Arctic/Siberian snowfall associated with sea-ice retreat [106] should be considered. As global warming continues, changes in snowfall patterns and their influence on the snow-AO teleconnection may become more pronounced [9], potentially enhancing heavy snowfall in the Alps. This trend is noticeable under emission scenarios RCP 2.6, 4.5 and 8.5, as indicated by Euro-CORDEX projections [107]. These observations align with an increasing trend of significant snowfall events across the Northern Hemisphere [3].

## 4.2. Model performance, limitations and forecasting errors

A key element in our forecast is the CFI, which combines AO, DMI, and TDsw and was used as an exogenous variable in both the location and log-scale components – Eqs. (2) and (3). However, the future values of CFI depend on projections of its individual components. These were forecasted using a PARMA(2,1) model, assuming periodic behaviour with a 25-year cycle. While this structure performed well under historical validation and passed heteroskedasticity testing (p=0.911), its application to future projections introduces a major limitation. As emphasised by Fereday et al. [14], the reliability of simulations of teleconnection patterns like AO and DMI is limited, especially under the influence of an accelerating Arctic amplification and shifting SST gradients. These indices are influenced by both internal variability and long-term forcing, which may lead to non-stationary behaviour that violates core assumptions of autoregressive models. Consequently, our model – while effective in capturing observed dynamics – shares the same shortcomings identified in the introduction to the article. Namely, purely statistical models can fail to capture structural changes in climate drivers, which could undermine the validity of long-term NHSD forecasts. To mitigate these risks, we constructed the CFI as a composite indicator, integrating multiple drivers to reduce over-reliance on any single index. While the empirical design of the CFI – Eq. (4) – improves correlation with NHSD (r=0.67, p<0.05), it still inherits uncertainty from its inputs. Despite this, our framework incorporates multiple layers of uncertainty quantification, including bootstrapped intervals and standard deviation windows. Moreover, the model performed well during validation, capturing key trends despite a non-random training split.

Despite the projected negative impact of global warming on snowfall in Northern Hemisphere mid-latitudes, heavy snowfall continue to play a key role in extreme weather events [3]. Understanding time-varying autoregressive models with input climate data is crucial for improving our predictive capabilities regarding snowfall patterns. However, forecasting errors are expected to increase when comparing forecasts for successive years, as they may not account for
<figure><img src="image_17.png" /></figure>

multi-decadal oscillations or changes in input variables like hydrological data, which inherently carry climate memory (e.g., persistence of prior climate conditions). Current approaches, such as the PARMAX(TVAR)-CSD method, struggle to accurately capture the changing frequency of severe weather events [108,109]. Moreover, while our method uses observational records, it may not encompass all possible weather patterns, resulting in inaccuracies [3], particularly in regions (e.g., western and southern Europe), where increasing temperatures reduce snowfall likelihood [110] and are projected to do so in the future [111]. The validation results highlight both the strengths and limitations of the PARMAX(TVAR)-CSD model. While the model effectively captures NHSD values – both within the narrower range of about 10–30 days observed over the last three decades and the larger fluctuations seen in the full time-series since 1883 – its reliance on large-scale indices makes it less adept to resolve fine-scale features (where thermodynamic and topographical complexities are likely to play a more significant role).

Although uncertainties persist regarding the impact of snowfall changes across different altitudes and potential transition from snowfall to rainfall, our model provides a credible foundation for forecasting NHSD dynamics, supported by rigorous validation and robust error statistics. Regarding the data split, while a random selection for both training and validation datasets would ensure a more balanced representation of historical climate variability, including the recent emergence of the anthropogenic forcing signal [112], the limited availability of historical data necessitated this approach. Although a comparative analysis focusing on periods with more pronounced anthropogenic impacts was not feasible due to the limitations of the time-series, the primary aim of the study was to assess the utility and predictive power of indices like the AO and DMI, combined with localised temperature shifts. Random data selection was intentionally avoided. Nevertheless, the model demonstrated robust performance during the validation period, and effectively captured key trends and patterns in the subsequent forecast period. This highlights its reliability and potential for accurate forecasting, suggesting that the model can effectively learn from the available data and provide meaningful predictions despite the non-random data split.

Nevertheless, future work should explore hybrid statistical-physical approaches that incorporate teleconnection forecasts from global climate models (e.g., Euro-CORDEX AO/DMI analogues). This would allow the model to adapt dynamically to evolving climatic conditions and improve resilience against potential structural breaks in atmospheric behaviour. In addition, while our dataset constraints limited the ability to assess the full impact of anthropogenic forcing periods separately, the model's successful capture of recent trends supports its current utility. Still, the dependence on autoregressive techniques highlights the importance of ongoing recalibration and validation as more data become available and as the climate system continues to evolve.

### 4.3. Socioeconomic implications of forecasted snowfall patterns

Model results suggest a strong temporal dependence in NHSD, shaped by both internal autocorrelation structures and external climate factors captured through the Climate Forcing Index (CFI). As the CFI reflects polar pressure variations, oceanic dipole patterns and seasonal temperature contrasts, its influence likely extends beyond the SPAR to the wider Swiss and Alpine regions. However, the spatial extent of CFI-driven changes and their broader effects on the Central European winter hydroclimate remain areas for future research. Forecasted declines in NHSD until around 2040 – supported by a statistically non-significant upward trend post-2035 – carry clear implications for winter tourism and hydro-power sectors, key economic drivers in the Alps that largely depend on snow cover. In Switzerland alone, passenger transport revenues during the winter season averaged €700 million between 2009 and 2013 [113]. Although this study did not assess sub-annual NHSD variability, the observed interannual fluctuations and downward trend suggest increasingly risk to critical tourism windows, such as the key Christmas/New Year holiday period. These changes also overlap with the annual low in natural river flow, posing risks to hydropower generation, which historically relied on snowmelt-fed runoff [113]. While warmer conditions may increase winter runoff from Alpine catchments [114], potentially offsetting some hydro-power concerns, shifts in the timing and form of precipitation could affect water management and downstream ecosystems [115]. Although the snow-to-rain transition was not directly modelled, reduced snow accumulation could still impact
groundwater recharge and long-term water storage, particularly in glacier-fed basins [71]. Looking ahead, projected NHSD recovery post-2045 – especially in higher elevations – may support a rebound in winter tourism and water availability, assuming other climatic and socioeconomic conditions remain stable.

### 4.4. Prospects for future research

The integration of autocorrelation function (ACF), partial autocorrelation function (PACF), Hurst exponent and attractor dynamics provides a robust analytical framework for understanding NHSD dynamics and addressing key challenges. While the model was not directly compared to one excluding these methods, their primary role is to inform model selection by identifying dependencies, long-term memory ($H=0.77$) and structural stability. These insights support models that incorporate both periodic and exogenous influences, enhancing interpretability rather than focusing solely on reducing predictive error. Collectively, these methodologies enhance predictive models by identifying dependencies, quantifying long-term memory and characterising the resilience and predictability of the system. In particular, ACF and PACF plots identify key time-series dependencies, informing model structure. The Hurst exponent analysis underscores trend persistence, emphasising long-term memory. Attractor analysis further highlights the system's structured and predictable behaviour, suggesting resilience to disturbances. These insights validate the use of the PARMAX(TVAR)-CSD model, which integrates periodicity, external predictors and memory, to provide accurate forecasts. Nonetheless, limitations remain, particularly in addressing novel climate scenarios and long-term shifts in extreme weather patterns.

One major limitation is that anthropogenic climate forcings such as greenhouse gas emissions and anthropogenic aerosols were not explicitly included in the analysis. While their effects may have been indirectly captured through historical climate patterns, especially with the extended training period covering 1884–2023, the direct causal relationships between anthropogenic forcings and snowfall are likely underrepresented. Given the complex and non-linear nature of these interactions [116], the model may not fully account for their influence, especially under extreme conditions. In addition, the uncertainty trajectory of anthropogenic forcings [117] complicates predictions for far-future scenarios (beyond 2060).

While the prediction of snowfall days in this study is based on historical climate data and atmospheric indicators, the incorporation of global climate models and ensemble techniques could improve the accuracy of forecasts and broaden their applicability across regions. However, caution should be exercised when extending hydrological predictions to longer time horizons due to significant uncertainties surrounding future precipitation estimates. Future research could also investigate underexplored or unidentified atmospheric modes of variability, and further refine existing models to improve predictive accuracy. However, the possibility to generalise the present study remains limited, as the model was validated solely with SPAR data. Expanding the analysis to include independent datasets from neighbouring alpine regions, such as the French and Italian Alps, could enhance the robustness of findings. A broader evaluation across diverse geographical contexts would provide valuable insights into the model's applicability beyond SPAR and inform adaptive strategies for managing snowfall-related impacts across the Alps. In addition, machine learning approaches may be worth exploring in future studies to identify non-linear patterns in the data that traditional methods might miss, further improving the model's accuracy and adaptability to complex climate change scenarios.

### Supporting information

**S1 Text. Climate Forcing Index forecasting.**
(DOCX)

**S1 Table. Complete set of observed and forecasted snowfall data that underpin the conclusions of this study.**
(XLSX)
S1 Fig. Climate Forcing Index (CFI) data (1884–2022). a. Annual CFI time-series with a long-term linear trend; b. Residuals of the CFI; c. Variance estimates of CFI with an 11-year moving window with linear trend regression line.

(TIF)

S2 Fig. a. Modelled (orange line) and observed (blue line) Climate Forcing Index (CFI) time-series during the validation stage (1994–2023); b. Observed (blue line) and modelled (orange line) CFI time-series during the forecasting stage (2024–2060). Training periods: a. 1884–1993; b. 1884–2023.

(TIF)

## Acknowledgments

We thank Folke Lundström (Stockholm) for help with formatting the references.

## Author contributions

**Conceptualization:** Nazzareno Diodato, Gianni Bellocchi.

**Data curation:** Nazzareno Diodato, Gianni Bellocchi.

**Formal analysis:** Nazzareno Diodato, Gianni Bellocchi.

**Funding acquisition:** Fredrik Charpentier Ljungqvist.

**Investigation:** Nazzareno Diodato, Gianni Bellocchi.

**Methodology:** Nazzareno Diodato, Gianni Bellocchi.

**Project administration:** Nazzareno Diodato, Gianni Bellocchi.

**Software:** Nazzareno Diodato, Gianni Bellocchi.

**Validation:** Nazzareno Diodato, Gianni Bellocchi.

**Visualization:** Nazzareno Diodato, Gianni Bellocchi.

**Writing – original draft:** Nazzareno Diodato, Fredrik Charpentier Ljungqvist, Gianni Bellocchi.

**Writing – review & editing:** Nazzareno Diodato, Fredrik Charpentier Ljungqvist, Gianni Bellocchi.

## References

1. Huss M, Bookhagen B, Huggel C, Jacobsen D, Bradley RS, Clague JJ, et al. Toward mountains without permanent snow and ice. Earth's Future. 2017;5(5):418–35. https://doi.org/10.1002/2016ef000514

2. McCrary RR, Mearns LO, Abel MR, Biner S, Bukovsky MS. Projections of North American snow from NA-CORDEX and their uncertainties, with a focus on model resolution. Clim Change. 2022;170(3–4):20. https://doi.org/10.1007/s10584-021-03294-8

3. Quante L, Willner SN, Middelanis R, Levermann A. Regions of intensification of extreme snowfall under future warming. Sci Rep. 2021;11(1):16621. https://doi.org/10.1038/s41598-021-95979-4 PMID: 34404852

4. Ombadi M, Risser MD, Rhoades AM, Varadharajan C. A warming-induced reduction in snow fraction amplifies rainfall extremes. Nature. 2023;619(7969):305–10. https://doi.org/10.1038/s41586-023-06092-7 PMID: 37380773

5. Preece JR, Mote TL, Cohen J, Wachowicz LJ, Knox JA, Tedesco M, et al. Summer atmospheric circulation over Greenland in response to Arctic amplification and diminished spring snow cover. Nat Commun. 2023;14(1):3759. https://doi.org/10.1038/s41467-023-39466-6 PMID: 37353537

6. Laurent L, Buoncristiani J-F, Pohl B, Zekollari H, Farinotti D, Huss M, et al. The impact of climate change and glacier mass loss on the hydrology in the Mont-Blanc massif. Sci Rep. 2020;10(1):10420. https://doi.org/10.1038/s41598-020-67379-7 PMID: 32591640

7. Notarnicola C. Overall negative trends for snow cover extent and duration in global mountain regions over 1982–2020. Sci Rep. 2022;12(1):13731. https://doi.org/10.1038/s41598-022-16743-w PMID: 35962171

8. Zhao K, Peng D, Gu Y, Luo X, Pang B, Zhu Z. Temperature lapse rate estimation and snowmelt runoff simulation in a high-altitude basin. Sci Rep. 2022;12(1):13638. https://doi.org/10.1038/s41598-022-18047-5 PMID: 35948622
9. Henderson GR, Peings Y, Furtado JC, Kushner PJ. Snow-atmosphere coupling in the Northern Hemisphere. Nature Clim Change. 2018;8(11):954–63. https://doi.org/10.1038/s41558-018-0295-6

10. Beniston M. Environmental Change in Mountains and Uplands. London: Arnold Publishers; 2000. 172 pp.

11. Tramblay Y, Koutroulis A, Samaniego L, Vicente-Serrano SM, Volaire F, Boone A, et al. Challenges for drought assessment in the Mediterranean region under future climate scenarios. Earth Sci Rev. 2020;210:103348. https://doi.org/10.1016/j.earscirev.2020.103348

12. Koehler J, Bauer A, Dietz AJ, Kuenzer C. Towards Forecasting Future Snow Cover Dynamics in the European Alps – The Potential of Long Optical Remote-Sensing Time Series. Remote Sensing. 2022;14(18):4461. https://doi.org/10.3390/rs14184461

13. Schiemann R, Athanasiadis P, Barriopedro D, Doblas-Reyes F, Lohmann K, Roberts MJ. Northern Hemisphere blocking simulation in current climate models: evaluating progress from the Climate Model Intercomparison Project Phase 5 to 6 and sensitivity to resolution. Weather Clim Dynam. 2020;1:277–92. https://doi.org/10.5194/wcd-1-277-202

14. Fereday D, Chadwick R, Knight J, Scaife AA. Atmospheric Dynamics is the Largest Source of Uncertainty in Future Winter European Rainfall. J Clim. 2018;31(3):963–77. https://doi.org/10.1175/jcli-d-17-0048.1

15. Puy A, Beneventano P, Levin SA, Lo Piano S, Portaluri T, Saltelli A. Models with higher effective dimensions tend to produce more uncertain estimates. Sci Adv. 2022;8(42):eabn9450. https://doi.org/10.1126/sciadv.abn9450 PMID: 36260678

16. Kirtman B, Power SB, Adedoyin JA, Boer GJ, Bojariu R, Camilloni J, et al. Near-term Climate Change: Projections and Predictability. In: 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, 2013.

17. Yoden S, Yoshida K. Impacts of Solar Activity Variations on Climate. In: Kusano K, editor. Solar-Terrestrial Environmental Prediction. Cham, Springer, 2023: 445–459. https://doi.org/10.1007/978-981-19-7765-7_15

18. Meredith M, Sommerkorn M, Cassotta S. Polar Regions. In: IPCC Special Report on the Ocean and Cryosphere in a Changing Climate. Cambridge, Cambridge University Press, 2019. https://doi.org/10.1017/9781009157964.005

19. Flanner MG, Shell KM, Barlage M, Perovich DK, Tschudi MA. Radiative forcing and albedo feedback from the Northern Hemisphere cryosphere between 1979 and 2008. Nature Geosci. 2011;4(3):151–5. https://doi.org/10.1038/ngeo1062

20. Deser C, Hurrell JW, Phillips AS. The role of the North Atlantic Oscillation in European climate projections. Clim Dyn. 2017;49(9–10):3141–57. https://doi.org/10.1007/s00382-016-3502-z

21. Marty C. Regime shift of snow days in Switzerland. Geophysical Research Letters. 2008;35(12):L12501. https://doi.org/10.1029/2008gl033998

22. Durand Y, Giraud G, Laternser M, Etchevers P, Mérindol L, Lesaffre B. Reanalysis of 47 Years of Climate in the French Alps (1958–2005): Climatology and Trends for Snow Cover. J Appl Meteorol Climatol. 2009;48(12):2487–512. https://doi.org/10.1175/2009jamc1810.1

23. Scherrer SC, Wüthrich C, Croci-Maspoli M, Weingartner R, Appenzeller C. Snow variability in the Swiss Alps 1864–2009. Intl Journal of Climatology. 2013;33(15):3162–73. https://doi.org/10.1002/joc.3653

24. Beniston M. Variations of snow depth and duration in the Swiss Alps over the last 50 years: links to changes in large-scale climatic forcings. Clim Chang. 1997;36:281–300. https://doi.org/10.1023/A:1005310214361

25. Irannezhad M, Ronkanen A-K, Kløve B. Wintertime climate factors controlling snow resource decline in Finland. Int J Climatol. 2016;36(1):110–31. https://doi.org/10.1002/joc.4332

26. Diodato N, Bellocchi G. Climate control on snowfall days in peninsular Italy. Theor Appl Climatol. 2020;140(3–4):951–61. https://doi.org/10.1007/s00704-020-03136-0

27. Bonsoms J, Salvador-Franch F, Oliva M. Snowfall and snow cover evolution in the Eastern Pre-Pyrenees (NE Iberian Peninsula). Geophys Res Lett. 2021;47(2):291–307. https://doi.org/10.18172/cig.4879

28. Annella C, Budillon G, Capozzi V. On the role of local and large-scale atmospheric variability in snow cover duration: a case study of Montevergine Observatory (Southern Italy). Environ Res Commun. 2023;5(3):031005. https://doi.org/10.1088/2515-7620/acc3e3

29. Sandgathe S, Brown BR, Carman JC, Infanti JM, Johnson B, McCarren D, et al. Exploring the Need for Reliable Decadal Prediction. Bulletin of the American Meteorological Society. 2020;101(2):E141–5. https://doi.org/10.1175/bams-d-19-0248.1

30. Gutowski WJ, Ullrich PA, Hall A, Leung LR, O'Brien TA, Patricola CM, et al. The Ongoing Need for High-Resolution Regional Climate Models: Process Understanding and Stakeholder Information. Bulletin of the American Meteorological Society. 2020;101(5):E664–83. https://doi.org/10.1175/bams-d-19-0113.1

31. Solari G. Le bucoliche di Publio Virgilio Marone recate in altrettanti versi italiani. Genova: Dalla Stamperia di Giossi; 1810. 288 pp.

32. Crumley CL. Historical Ecology: A Robust Bridge between Archaeology and Ecology. Sustainability. 2021;13(15):8210. https://doi.org/10.3390/su13158210

33. Chatfield C. Time-series forecasting. 1st edition. New York: Chapman and Hall/CRC; 2000.

34. Mulligan M, Wainwright J. Modelling and model building. In: Wainwright J, Mulligan M, editors. Environmental modelling: finding simplicity in complexity. New York: John Wiley & Sons, Ltd; 2004. pp. 7–73.

35. Cohen J, Rind D. The effect of snow cover on the climate. J Clim. 1991;4:689–706. https://doi.org/10.1175/1520-0442(1991)004<0689:TEOSCO≥2.0.CO;2
36. Dutra E, Schär C, Viterbo P, Miranda PMA. Land–atmosphere coupling associated with snow cover. Geophysical Research Letters. 2011;38(15):L15707. https://doi.org/10.1029/2011gl048435

37. Derksen C, LeDrew E, Goodison B. SSM/I derived snow water equivalent data: The potential for investigating linkages between snow cover and atmospheric circulation. Atmosphere-Ocean. 2010;36(2):95–117. https://doi.org/10.1080/07055900.1998.9649608

38. Carleton TA, Hsiang SM. Social and economic impacts of climate. Science. 2016;353(6304):aad9837. https://doi.org/10.1126/science.aad9837 PMID: 27609899

39. Ribes A, Qasmi S, Gillett NP. Making climate projections conditional on historical observations. Sci Adv. 2021;7(4):eabc0671. https://doi.org/10.1126/sciadv.abc0671 PMID: 33523939

40. Zhang Y, Yang H, Cui H, Chen Q. Comparison of the Ability of ARIMA, WNN and SVM Models for Drought Forecasting in the Sanjiang Plain, China. Nat Resour Res. 2020;29(2):1447–64. https://doi.org/10.1007/s11053-019-09512-6

41. Islam F, Imteaz MA. Use of Teleconnections to Predict Western Australian Seasonal Rainfall Using ARIMAX Model. Hydrology. 2020;7(3):52. https://doi.org/10.3390/hydrology7030052

42. Ponnamperuma N, Rajapakse L. Comparison of Time Series Forecast Models for Rainfall and Drought Prediction. 2021 Moratuwa Engineering Research Conference (MERCon) 27–29 July 2021, IEEE Xplore Moratuwa, Sri Lanka. Available from: https://doi.org/10.1109/MERCon52712.2021.9525690

43. Meng L, Zhu L. Statistical Modeling of Monthly and Seasonal Michigan Snowfall Based on Machine Learning: A Multiscale Approach. Artificial Intelligence for the Earth Systems. 2022;2(4):e230016. https://doi.org/10.1175/aies-d-23-0016.1

44. Monteiro D, Morin S. Multi-decadal analysis of past winter temperature, precipitation and snow cover data in the European Alps from reanalyses, climate models and observational datasets. Cryosphere. 2023;17(8):3617–60. https://doi.org/10.5194/tc-17-3617-2023

45. Clark CA, Elless TJ, Lyza AW, Ganesh-Babu B, Koning DM, Carne AR, et al. Spatiotemporal Snowfall Variability in the Lake Michigan Region: How is Warming Affecting Wintertime Snowfall? J Appl Meteorol Climatol. 2016;55(8):1813–30. https://doi.org/10.1175/jamc-d-15-0285.1

46. McCray CD, Schmidt GA, Paquin D, Leduc M, Bi Z, Radiyat M, et al. Changing Nature of High-Impact Snowfall Events in Eastern North America. J Geophys Res. 2023;128(13):e2023JD038804. https://doi.org/10.1029/2023jd038804

47. Hammer HL. Statistical models for short- and long-term forecasts of snow depth. Journal of Applied Statistics. 2017;45(6):1133–56. https://doi.org/10.1080/02664763.2017.1357683

48. MeteoSwiss. Federal Office of Meteorology and Climatology. 2024. Available from: https://www.meteoswiss.admin.ch/climate/climate-change/heat-waves-droughts-cold-and-snowfall.html

49. Goudie AS. Global warming and fluvial geomorphology. Geomorphology. 2006;79(3–4):384–94. https://doi.org/10.1016/j geomorph.2006.06.023

50. Beniston M, Goyette S. Changes in variability and persistence of climate in Switzerland: Exploring 20th century observations and 21st century simulations. Global and Planetary Change. 2007;57(1–2):1–20. https://doi.org/10.1016/j.gloplacha.2006.11.004

51. Cache T, Ramirez JA, Molnar P, Ruiz-Villanueva V, Peleg N. Increased erosion in a pre-Alpine region contrasts with a future decrease in precipitation and snowmelt. Geomorphology. 2023;436:108782. https://doi.org/10.1016/j geomorph.2023.108782

52. Diodato N, Fratianni S, Bellocchi G. Reconstruction of snow days based on monthly climate indicators in the Swiss pre-alpine region. Reg Environ Change. 2020;20(2):55. https://doi.org/10.1007/s10113-020-01639-0

53. Lit R. Forecasting the VIX in the midst of COVID-19. Time Series Lab – Article Series, 2020. 15 pp. Available from: https://timeserieslab.com/articles/rlit_vix_v2.pdf

54. Schär S, Marelli S, Sudret B. Emulating the dynamics of complex systems using autoregressive models on manifolds (mNARX). Mech Syst Signal Process. 2024;208:110956. https://doi.org/10.1016/j.ymssp.2023.110956

55. Antoine E, Marquer L, Muigg B, Tegel W, Bisson U, Bolliger M, et al. Legacy of last millennium timber use on plant cover in Central Europe: Insights from tree rings and pollen. Sci Total Environ. 2024;922:171157. https://doi.org/10.1016/j.scitotenv.2024.171157 PMID: 38412879

56. Joly D, Castel T, Pohl B, Richard Y. Influence of spatial information resolution on the relation between elevation and temperature. Intl Journal of Climatology. 2018;38(15):5677–88. https://doi.org/10.1002/joc.5771

57. Faranda D. An attempt to explain recent changes in European snowfall extremes. Weather Clim Dynam. 2020;1(2):445–58. https://doi.org/10.5194/wcd-1-445-2020

58. Cattiaux J, Vautard R, Cassou C, Yiou P, Masson-Delmotte V, Codron F. Winter 2010 in Europe: A cold extreme in a warming climate. Geophys Res Lett. 2010;37(20):L20704. https://doi.org/10.1029/2010gl044613

59. Hartwig F, Dearing BE. Exploratory data analysis. Newbury Park: SAGE University Paper; 1979, 83 pp.

60. Franzke CLE, Barbosa S, Blender R, Fredriksen H, Laepple T, Lambert F, et al. The Structure of Climate Variability Across Scales. Reviews of Geophysics. 2020;58(2):e2019RG000657. https://doi.org/10.1029/2019rg000657

61. Ghil M, Allen MR, Dettinger MD, Ide K, Kondrashov D, Mann ME, et al. Advanced spectral methods for climatic time series. Rev Geophys. 2002;40(1):3.1–3.41. https://doi.org/10.1029/2000rg000092

62. Krawiecki A, Holyst JA, Helbing D. Volatility clustering and scaling for financial time series due to attractor bubbling. Phys Rev Lett. 2002;89(15):158701. https://doi.org/10.1103/PhysRevLett.89.158701 PMID: 12366029
63. Kirichenko L, Radivilova T, Bulakh V. Generalized approach to hurst exponent estimating by time series. Informatyka, Automatyka, Pomiary W Gospodarce I Ochronie Środowiska. 2018;8(1):28–31. https://doi.org/10.5604/01.3001.0010.8639

64. Ramachandra AR, Bhattacharya D. Effect of short-term memory on Hurst phenomenon. J Hydrol Eng. 2001;6:125–31. https://doi.org/10.1061/(ASCE)1084-0699(2001)6:2(125)

65. Valenzuela O, Rojas I, Rojas F, Pomares H, Herrera LJ, Guillen A, et al. Hybridization of intelligent techniques and ARIMA models for time series prediction. Fuzzy Sets and Systems. 2008;159(7):821–45. https://doi.org/10.1016/j.fss.2007.11.003

66. McLeod AI, Li WK. Diagnostic checking for arma time series models using squared-residual autocorrelations. J Time Ser Anal. 1983;4(4):269–73. https://doi.org/10.1111/j.1467-9892.1983.tb00373.x

67. Mandelbrot BB, Taqqu MS. Robust R/S analysis of long run serial correlation. Bull Int Stat Inst. 1979;48:59–104.

68. Taheriyoun AR, Moghimbeygi M. Visual information and expert's idea in Hurst index estimation of the fractional Brownian motion using a diffusion type approximation. Sci Rep. 2017;7:42482. https://doi.org/10.1038/srep42482 PMID: 28195153

69. Roman HE, Porto M. Fractional derivatives of random walks: time series with long-time memory. Phys Rev E Stat Nonlin Soft Matter Phys. 2008;78(3 Pt 1):031127. https://doi.org/10.1103/PhysRevE.78.031127 PMID: 18851013

70. Montgomery DC, Jennings CL, Kulahci M. Introduction to time series analysis and forecasting. John Wiley & Sons. 2016. 627 p.

71. Li Z, Chen Y, Li Y, Wang Y. Declining snowfall fraction in the alpine regions, Central Asia. Sci Rep. 2020;10(1):3476. https://doi.org/10.1038/s41598-020-60303-z PMID: 32103068

72. Koopman SJ, Lucas A, Scharth M. Predicting Time-Varying Parameters with Parameter-Driven and Observation-Driven Models. Rev Econ Stat. 2016;98(1):97–110. https://doi.org/10.1162/rest_a_00533

73. Basawa IV, Lund R, Shao Q. First-order seasonal autoregressive processes with periodically varying parameters. Stat Probab Lett. 2004;67(4):299–306. https://doi.org/10.1016/j.spl.2004.02.001

74. Lit R, Koopman SJ, Harvey AC. Time Series Lab, Score Edition. 2022.

75. Sun Y, Wang S-YS, Li R, Buckley BM, Gillies R, Hansen KG. Feasibility of Predicting Vietnam's Autumn Rainfall Regime Based on the Tree-Ring Record and Decadal Variability. Climate. 2018;6(2):42. https://doi.org/10.3390/cli6020042

76. Liu Y, Donat MG, England MH, Alexander LV, Hirsch AL, Delgado-Torres C. Enhanced multi-year predictability after El Niño and La Niña events. Nat Commun. 2023;14(1):6387. https://doi.org/10.1038/s41467-023-42113-9 PMID: 37821438

77. Lüthi S, Ban N, Kotlarski S, Steger CR, Jonas T, Schär C. Projections of Alpine Snow-Cover in a High-Resolution Climate Simulation. Atmosphere. 2019;10(8):463. https://doi.org/10.3390/atmos10080463

78. Bartolini WM. Convection-permitting ensemble forecasts of the 10–12 December 2013 lake-effect snow event: Sensitivity to microphysical, planetary boundary layer, and surface layer parameterizations. 2019. Ph.D. thesis, State University of New York at Albany. Available from: https://doi.org/10.54014/JV49-JRFZ

79. Niziol TA, Snyder WR, Waldstreicher JS. Winter Weather Forecasting throughout the Eastern United States. Part IV: Lake Effect Snow. Wea Forecasting. 1995;10(1):61–77. https://doi.org/10.1175/1520-0434(1995)010<0061:wwfte>2.0.co;2

80. Zhang T, Zhu S, Song Y, Wang X, Chen H. The Leading Modes of Northern Eurasian Winter Snowfall Variability and the Potential Influencing Factors. Journal of Climate. 2023;36(22):7811–26. https://doi.org/10.1175/jcli-d-22-0881.1

81. Liu Y, Wang L, Zhou W, Chen W. Three Eurasian teleconnection patterns: spatial structures, temporal variability, and associated winter climate anomalies. Clim Dyn. 2014;42(11–12):2817–39. https://doi.org/10.1007/s00382-014-2163-z

82. Jiang Y, Cheung H-N, Li Y, Yang S. Intra-seasonal variation of the wintertime Polar/Eurasia pattern. Clim Dyn. 2022;61(1–2):813–30. https://doi.org/10.1007/s00382-022-06612-7

83. Haben S. Verification and Evaluation of Load Forecast Models. In: Haben S, Voss M, Holderbaum W, editors. Core Concepts and Methods in Load Forecasting, Springer Cham; 2023, p. 89–105. https://doi.org/10.1007/978-3-031-27852-5_7

84. Trouet V, Van Oldenborgh GJ. KNMI Climate Explorer: A Web-Based Research Tool for High-Resolution Paleoclimatology. Tree-Ring Research. 2013;69(1):3–13. https://doi.org/10.3959/1536-1098-69.1.3

85. Bartolini E, Claps P, D'Odorico P. Connecting European snow cover variability with large scale atmospheric patterns. Adv Geosci. 2010;26:93–7. https://doi.org/10.5194/adgeo-26-93-2010

86. Yuan C, Tozuka T, Miyasaka T, Yamagata T. Respective influences of IOD and ENSO on the Tibetan snow cover in early winter. Clim Dyn. 2009;33(4):509–20. https://doi.org/10.1007/s00382-008-0495-2

87. Scherrer S, Appenzeller C. Swiss Alpine snow pack variability: major patterns and links to local climate and large-scale flow. Clim Res. 2006;32:187–99. https://doi.org/10.3354/cr032187

88. Pettitt AN. A Non-Parametric Approach to the Change-Point Problem. Appl Stat. 1979;28(2):126. https://doi.org/10.2307/2346729

89. North RP, Livingstone DM, Hari RE, Köster O, Niederhauser P, Kipfer R. The physical impact of the late 1980s climate regime shift on Swiss rivers and lakes. Inland Waters 2013;3:341–350. https://doi.org/10.5268/IW-3.3.560

90. Reid PC, Hari RE, Beaugrand G, Livingstone DM, Marty C, Straile D, et al. Global impacts of the 1980s regime shift. Glob Chang Biol. 2016;22(2):682–703. https://doi.org/10.1111/gcb.13106 PMID: 26598217
91. Shrestha DL, Solomatine DP. Data-driven approaches for estimating uncertainty in rainfall-runoff modelling. International Journal of River Basin Management. 2008;6(2):109-22. https://doi.org/10.1080/15715124.2008.9635341

92. Lim KJ, Engel BA, Tang Z, Muthukrishnan S, Choi J, Kim K. Effects of calibration on L-THIA GIS runoff and pollutant estimation. J Environ Manage. 2006;78(1):35-43. https://doi.org/10.1016/j.jenvman.2005.03.014 PMID: 16112801

93. Diodato N, Bellocchi G, Bertolin C, Camuffo D. Mixed nonlinear regression for modelling historical temperatures in Central-Southern Italy. Theor Appl Climatol. 2013;113(1-2):187-96. https://doi.org/10.1007/s00704-012-0775-y

94. Diodato N, Mao L, Borrelli P, Panagos P, Fiorillo F, Bellocchi G. Climate-scale modelling of suspended sediment load in an Alpine catchment debris flow (Rio Cordon-northeastern Italy). Geomorphology. 2018;309:20-8. https://doi.org/10.1016/j geomorph.2018.02.026

95. Box GEP, Jenkins G. Time series analysis: Forecasting and control. San Francisco: Holden-Day; 1970.

96. Sun Y, Solomon S, Dai A, Portmann RW. How Often Does It Rain?. Journal of Climate. 2006;19(6):916-34. https://doi.org/10.1175/jcli3672.1

97. Hyndman RJ, Athanasopoulos G. Forecasting: principles and practice. 3rd edition. Melbourne, Australia: OTexts; 2020. Available from: https://otexts.com/fpp2

98. Swanson DA, Tayman J, Bryan TM. MAPE-R: a rescaled measure of accuracy for cross-sectional subnational population forecasts. J Pop Research. 2011;28(2-3):225-43. https://doi.org/10.1007/s12546-011-9054-5

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

100. IPCC. Synthesis Report of the Sixth Assessment Report (AR6) Longer Report; Intergovernmental Panel on Climate Change: Geneva, Switzerland, 2022. Available from: https://report.ipcc.ch/ar6syr/pdf/IPCC_AR6_SYR_LongerReport.pdf

101. Burgess MG, Ritchie J, Shapland J, Pielke R Jr. IPCC baseline scenarios have over-projected CO₂ emissions and economic growth. Environ Res Lett. 2020;16(1):014016. https://doi.org/10.1088/1748-9326/abcdd2

102. Cao W, Duan C, Yang T, Wang S. Higher Heat Stress Increases the Negative Impact on Rice Production in South China: A New Perspective on Agricultural Weather Index Insurance. Atmosphere. 2022;13(11):1768. https://doi.org/10.3390/atmos13111768

103. Le Roux E, Evin G, Samacoïts R, Eckert N, Blanchet J, Morin S. Projection of snowfall extremes in the French Alps as a function of elevation and global warming level. The Cryosphere. 2023;17(11):4691-704. https://doi.org/10.5194/tc-17-4691-2023

104. Hind A, Zhang Q, Brattström G. Problems encountered when defining Arctic amplification as a ratio. Sci Rep. 2016;6:30469. https://doi.org/10.1038/srep30469 PMID: 27461918

105. Cohen J, Agel L, Barlow M, Entekhabi D. No detectable trend in mid-latitude cold extremes during the recent period of Arctic amplification. Commun Earth Environ. 2023;4(1):341. https://doi.org/10.1038/s43247-023-01008-9

106. Ghatak D, Frei A, Gong G, Stroeve J, Robinson D. On the emergence of an Arctic amplification signal in terrestrial Arctic snow extent. J Geophys Res. 2010;115(D24):D24105. https://doi.org/10.1029/2010jd014007

107. Larsson Ivanov O, Bärring L, Wilcke RAI. Climate change impact on snow loads in northern Europe. Structural Safety. 2022;97:02231. https://doi.org/10.1016/j.strusafe.2022.102231

108. Bavay M, Grünewald T, Lehning M. Response of snow cover and runoff to climate change in high Alpine catchments of Eastern Switzerland. Advances in Water Resources. 2013;55:4-16. https://doi.org/10.1016/j.advwatres.2012.12.009

109. Bard A, Renard B, Lang M, Giuntoli I, Korck J, Koboltschnig G, et al. Trends in the hydrologic regime of Alpine rivers. J Hydrol. 2015;529:1823-37. https://doi.org/10.1016/j.jhydrol.2015.07.052

110. Deser C, Phillips AS, Alexander MA, Smoliak BV. Projecting North American Climate over the Next 50 Years: Uncertainty due to Internal Variability*. J Clim. 2014;27(6):2271-96. https://doi.org/10.1175/jcli-d-13-00451.1

111. Jones AD, Rastogi D, Vahmani P, Stansfield AM, Reed KA, Thurber T, et al. Continental United States climate projections based on thermodynamic modification of historical weather. Sci Data. 2023;10(1):664. https://doi.org/10.1038/s41597-023-02485-5 PMID: 37770463

112. Risser MD, Collins WD, Wehner MF, O'Brien TA, Huang H, Ullrich PA. Anthropogenic aerosols mask increases in US rainfall by greenhouse gases. Nat Commun. 2024;15(1):1318. https://doi.org/10.1038/s41467-024-45504-8 PMID: 38388495

113. Schmucki E, Marty C, Fierz C, Weingartner R, Lehning M. Impact of climate change in Switzerland on socioeconomic snow indices. Theor Appl Climatol. 2017;127(3-4):875-89. https://doi.org/10.1007/s00704-015-1676-7

114. Fischer EM, Knutti R. Observed heavy precipitation increase confirms theory and early models. Nature Clim Change. 2016;6(11):986-91. https://doi.org/10.1038/nclimate3110

115. Prein AF, Heymsfield AJ. Increased melting level height impacts surface precipitation phase and intensity. Nat Clim Chang. 2020;10(8):771-6. https://doi.org/10.1038/s41558-020-0825-x

116. Chen H, Sun J, Lin W. Anthropogenic influence would increase intense snowfall events over parts of the Northern Hemisphere in the future. Environment Res Lett. 2020;15(11):114022. https://doi.org/10.1088/1748-9326/abbc93

117. Krasting JP, Broccoli AJ, Dixon KW, Lanzante JR. Future Changes in Northern Hemisphere Snowfall. Journal of Climate. 2013;26(20):7813-28. https://doi.org/10.1175/jcli-d-12-00832.1