1 Introduction

Quando Orion dal cielo

Declinando imperversa,

E pioggia e nevi e gelo

Sopra la terra ottenebrata versa,

me spinto ne la iniqua

stagione, infermo il piede,

tra il fango e tra l’obliqua

furia de’ carri la città gir vede

When Orion from heaven

Declining raging,

And rain and snow and frost

Over the darkened earth pours,

pushed me into the iniquitous

season, foot infirm,

between the mud and the oblique

fury of the chariots the town sees

Giuseppe Parini, La caduta, 1785, published by Reina 1825, p. 207

(our translation into English)

Injured by a fall on a snowy and icy ground, the Italian poet Giuseppe Parini (1729–1799) begins his ode La caduta (“The fall”) with the epic description of an eighteenth century winter season in his city (Milan, northern Italy). The Italians of the eighteenth century had certainly to face a much more difficult challenge than today when it came to facing winter. According to Stoppani (1876), between 1797 and 1806 the total number of days with snow in Milan had been 243 (i.e., 26 on average per year), but the situation had already substantially changed from 1857 to 1876 with a total of 166 snow days (i.e., eight snow days per year). Extensive evidence reveals in fact that Earth’s snow cover and snow amount have been decreasing, especially in the Northern Hemisphere, where this declining trend is expected to affect climate variability by the albedo effect (Beniston 1997, 2003; Mote et al. 2018; Brown et al. 2017; Fontrodona Bach et al. 2018; Pulliainen et al. 2020; Diodato et al. 2021). Although the importance of snowfall and snow cover on climate had already been highlighted in the nineteenth century (Woeikof 1885), it is only recently that the associated impacts on human livelihood, economy and ecosystems have been considered within the context of contemporary global warming (e.g., Xu et al. 2009; Beniston 2012) or studied for historical times (e.g., Ljungqvist et al. 2021). In particular, snow cover changes alter the diabatic heating of the Earth’s surface by increasing the fraction of solar radiation reflected away by the surface (i.e., the surface albedo), thus becoming an important constituent of the terrestrial radiation balance (e.g., Goward 2005; Déry and Brown 2007; Sandells and Flocco 2014). In particular, they are an important positive feedback of temperature amplification in the Arctic and snow-dominated high-elevation regions (Pepin et al. 2015). Snow cover also affects the timing and flow peaks produced by snow melting in catchment areas (Wang et al. 2015) by delaying the transfer of precipitation to surface runoff and infiltration, which influences flow reduction and the extent of the flow network during summer base flow (Yarnell et al. 2010; Godsey et al. 2014). Alterations in snow cover may, in turn, influence mountain ecosystem hydrology (Gobiet et al. 2014), biological and chemical processes in soil and water (Magnani et al. 2017), crop phenology and plant composition (Rogora et al. 2018). In this context, the inter-annual variability of cool-season precipitation is particularly important at middle to high latitudes, where winter precipitation likely occurs as snow at high elevations through orographic enhancement (e.g., Dettinger et al. 1998; Selkowitz et al. 2002; Terzago et al. 2012; Luce et al. 2013). Thus, the past frequency of snowfall occurrence and spatiotemporal trends of snow depth and snow cover extent are important determinants of the climate signal towards detecting and discovering potential changes of climate that may lead to considerable impacts on both the environment and society (O’Gorman 2014; Viru and Jaagus 2020). These determinants are also important for the validation of climate models and a better understanding of the interactions among the different spheres (geosphere, biosphere, atmosphere, hydrosphere and cryosphere) of the Earth’s system (Brown 2014).

Snow cover duration, i.e., number of days with snow cover on the ground, is the most suitable climatic characteristic to describe snow cover conditions. Snow cover data are currently gathered through ground-based observations and by satellite observations of the Earth’s surface (Robinson 1989; Estilow et al. 2015). However, observational datasets are still limited worldwide, with measurements of snow cover available from the late nineteenth century or the twentieth century in parts of the Northern Hemisphere (Brown 2000). For instance, snow cover duration observations started in the 1890s in Estonia (Jaagus 1997), in the 1920s in northeastern USA (Leathers and Luff, 1997), in the 1930s in Bulgaria (Brown and Petkova 2007), in the 1950s in the French Alps (Durand et al. 2009) and in western China (Dahe 2006), in the 1960s in the Qinghai-Tibetan Plateau (Xu 2017), and in the 1970s in the Swiss Alps (Klein et al. 2016). Snow cover records through satellite observations point to a broad pattern of decreasing snow cover, particularly in spring (Kunkel et al. 2016). The contribution of patchy (fractional) snow cover can be prominent in complex terrains but our ability to monitor trends in mountain regions appears limited. The coarse spatial resolution is still problematic for several mountain regions, with broad mountain systems only partially captured by these data, though higher-resolution new satellite missions are helping to fill this gap (Bormann et al. 2018). In Italy, published time-series of snow cover data extend back to the 1910s for the subalpine area (Mercalli and Cat Berro 2005), the 1930s for the Alps (western sector, Bovo et al. 2009; eastern sector, Cat Berro et al. 2008) and for the Po Valley (Parma observatory, Diodato et al. 2020a), and the 1950s for the northern Apennines (De Bellis et al. 2010).

For this reason, the study of past climate is important because it increases our understanding of the magnitude and patterns of long-term climate variability, including the direction and magnitude of changes of extremes such as excessive cold and snowy periods (Dong and Menzel 2020). In Europe, the snowfall dynamics are modulated by the latitude and atmospheric circulation patterns, and local orographic situations (Croce et al. 2018; Kretschmer et al. 2018). The Atlantic polar sea air masses collide here with the continental polar air masses associated with the Asian high pressure (Bednorz 2004). Temperature control of snow cover is dominant in snowy transition regions where the mean winter temperature is slightly lower than, and often crosses, the melting point (Diodato and Bellocchi 2020). The situation is different in regions where the mean winter temperature is considerably lower than 0 °C, for instance in the Alps, where changes in snow cover days are mainly controlled by precipitation regimes (Clark et al. 1999). The plains of northern Italy also experience a rather continental climate (Diodato et al. 2020a). Many of these plains were often covered by snow and ice in winter during past cold periods. At low European latitudes, snow cover duration is mainly affected by air temperature and only secondarily by precipitation conditions. However, at higher elevations, which are characterized by long periods with snow cover, precipitation plays a major role in driving snowfall (Diodato and Bellocchi 2020).

The presence of vast water reservoirs in southern Italy as a testimony of past colder climate conditions of the Little Ice Age (LIA, ~ 1300–1850 CE) is also an aspect highlighted in historical chronicles (Gentilcore 2019). For the cathedral canon Carlo Celano (1625–1693), water in Naples (40° 50′ N, 14° 15′ E) in the seventeenth century was so cold that it is difficult to believe that snow has not been added to it (Celano 1692, vol. 2, p. 203). The physician Salvatore De Renzi (1800–1872), in his report on the medical topography of the Kingdom of Naples in the nineteenth century (De Renzi 1829, p. 170), wrote “The snow falls in November above the inland mountains, and within the higher valley, and snow-cover lasts up to end of May.” Today’s emerging negative trends in snow depth and cover days suggest that snow-dependent regions in southern Europe could be affected by these changes because the number of snow days per year and the days with snow on the ground (DSG, day year−1) determine soil processes and water flow in rivers, streams, lakes and reservoirs (Diodato et al. 2019, 2020a, c, 2021). The Italian territory (the focus of this study, excluding the islands) is roughly delimited to the north by the alpine range (around 46° N) and its southern-most point is roughly 38° N, divided into east and west by the Apennine mountains. For this study, we made the hypothesis that DSG estimates across Italy can be obtained based on climate and elevation data from a limited set of sites or areas located around 41° N (i.e., in the region of the city of Naples). The 41st parallel North, which roughly marks the border of the prevalently Mediterranean climate characterising southern Italy, also crosses a sector of peninsular Italy where there is a rich heritage of snowfall data. This territory thus represents a strategic area for understanding and interpreting snowfall dynamics across the Italian peninsula (Diodato and Bellocchi 2020). This was done with the statistically based, physically consistent DSG (CSA) model (developed in this study), where CSA stands for central-southern Apennines, which include two sectors: the Samnite Apennines and the Abruzzi Apennines. The Samnite area identifies a subdivision of the southern Apennines, lying in the ancient region of Samnium, roughly between Molise and the northern sector of the Campania region. Contiguous to Samnium lies the Abruzzi Apennine chain, in the southern segments of the central Appenine mountains, which contains the most rugged terrain of the Apennines. The DSG model was developed within the CSA (first objective of the study) and used for an overview of scale concepts and to illustrate estimates of snow cover days of Italy (second objective of the study). In this way, we assessed the quality of DSG estimates and examined their spatio-temporal variation over complex terrains, within and beyond the CSA. A focus on the CSA was prepared by an analysis of the hydrology and climatology of the region.

2 Materials and methods

2.1 Environmental setting

With cool and wet winters, and hot and dry summers, the majority of Italy experiences a Mediterranean climate (Lionello et al. 2006). Orographically enhanced precipitation and highly variable temperatures determine the seasonality of the snow cover (Fayad et al. 2017). Under these conditions, the modelling of DSG is difficult. A parsimonious modelling approach was developed for the analysis of DSG patterns over two complex-terrain configurations of ~ 10,000 km2 each, placed at the core of the Mediterranean region (Fig. 1a, b). The Samnite area stretches across the Campania and Molise regions, around 41° 20′ N and 14° 25′ E, with elevations ranging from about 100 to 2050 m a.s.l. (Fig. 1c). The second area is located in the Abruzzi region, slightly further eastward (~ 13°E) and northward (above 42° N). The typical Mediterranean climate of this part of the Italian peninsula has more continental features than the maritime climate of the Tyrrhenian coast (e.g., the Naples area). This is mainly due to the topographical structure of central/southeastern Italy, where high-altitude mountain chains constitute a component of considerable extent.

Fig. 1
figure 1

a Geographical setting (orange square) and map of Italy, b calibration sites in the Samnite Apennines (red points) and Abruzzi Apennines areas (black square), and c validation sites in the Po River basin, and in the Samnite and Calabrian Apennines (red points), and areas (red rectangles) in western Alps and northern Apennines. BNV: Benevento Valley, MVOBS: Montevergine Observatory, MM: Matese Massif; ST: S. Tommaso OBS. Maps are authors’ own elaboration from free, public domain, images: a Wikimedia Commons (https://commons.wikimedia.org/wiki/File:BlankMap-Europe-v4.png); b and c Wikipedia (https://it.wikipedia.org/wiki/Italia)

Full size image

From the big mountains of the Matese Massif (MM) to the north (at the border between the Campania and Molise regions), through the Taburno Massif and the slopes of the Partenio mountains (which include Montevergine) to the south, the Samnium area manifests an intricate labyrinth of tangles and reliefs in a plurality of environments, settings and variations (Fig. 2). Mountains and hills of this internal landscape together cover about 80% of the area, which offers the panoramic view that was delineated in the description of De Renzi (1829, p. 12): “While the whole extension of the Mediterranean is in front of us and we are surrounded by a warm volcanic atmosphere, we are then surrounded to the south-east by a chain of the Apennines, we have at a straight distance of 20 miles high mountains like Taburno, and Partenio, and at no more than 30 miles radius, the chilling peaks of the Matese.”

Fig. 2
figure 2

Perspective view of Samnite Apennines between Montevergine and Benevento Valley (BNV) to the south, and Matese Massif to the north. Elevations in metres above sea level are reported for reference sites. Map is an output image created from OpenStreet Map (https://tinyurl.com/y4qa79fu)

Full size image

In both the massifs of Montevergine (up to some 1590 m a.s.l.) and Matese (up to some 2050 m a.s.l.), peaks are above the limit of the free troposphere (~ 1500 m a.s.l. or 850 hPa; Jones et al. 2003), with limited anthropogenic influences. In this inland landscape, the distribution of snow is extremely varied and irregular in relation to the orography and exposure of the slopes. In particular, average seasonal and annual snow depths, the frequency of the phenomenon, and the permanence of snow on the ground increase with the altitude and the slope of the site, as well as other factors such as the distance from the Tyrrhenian Sea. Climate in these ranges of the Italian hinterland is variously affected by two coastal areas, the Tyrrhenian to the west and the Adriatic to the east. The Adriatic side of the Samnium, toward the town of Bojano (41° 29′ N, 14° 28′ E, 482 m a.s.l.), is snowier, due to the limited mitigating effect of the Adriatic Sea, which being smaller and less deep is also colder than the Tyrrhenian Sea. Exposed to the north-east, this side of the Italian peninsula is under the influence of snow-bearing cold fronts from continental Europe. In the inner Abruzzi region, climate is markedly continental with typical mountain characteristics. Here, night-frosts are frequent, widespread and intense, with persistent snow on the ground and particularly intense cold waves, as highlighted in the description made by the Italian politician and scholar Giuseppe Del Re (1806–1864), e.g., Snows, which begin in October in the mountains, only liquefy towards the end of June (Del Re 1835, p. 14).

Synoptic situations favourable to snowfall occurrence and persistence in central/south-eastern Italy usually derive from a low-pressure system located over southern Italy and a high pressure centred over Scandinavia (Fazzini et al. 2005). Another configuration, less recurrent, is the occurrence of a low-pressure centre over the Ionian Sea (extreme south of Italy) extending from western Russia, passing over the Aegean Sea and influencing snowfall in the Mediterranean (Dafis et al. 2016). Drawing cold air from the Balkan Peninsula, such a low-pressure centre is a source of humidity, a necessary condition for snowfall. Snow episodes can also originate in situ after the advection of very cold air masses, which remain trapped in the Benevento Valley, above which mild air flows and produces snowfall, sometimes with great abundance. Figure 3 shows remarkably different seasonal and monthly patterns of snow depth between the Benevento Valley (BNV; Fig. 3a) and Montevergine Observatory (MVOBS, Fig. 3b).

Fig. 3
figure 3

Monthly mean distribution of snow depth (black curves) and seasonal precipitation regime (coloured bands with grey dotted curves) at a Benevento Valley (BNV), and b Montevergine Observatory (MVOBS), with related bounds of the 10th, 25th, 75th and 90th percentiles (dark and light grey envelops around the mean black curve), calculated over the period 1971–2018. The locations of the two sites are shown in Fig. 1b. The graphs are an output generated from data freely available in F4map Demo—Interactive 3D map (https://demo.f4map.com/#camera.theta=0.9)

Full size image

2.2 Snow cover data resources

A day with snow cover is a day with snow of at least 0.01 m depth observed in over 50% of an open neighbouring area (UNESCO/IASH/WMO 1970; WMO 2008). Winter maximum number of days with snow depth matching these criteria was used to characterize snow cover duration (e.g., Falarz 2004; Petkova et al. 2004). However, in the absence of actual observations for determining the correct values of snow cover duration in the past, the proofs of the existence of snow on the ground are those substantiated by documentary evidence. Studies aimed at collecting, reconstructing and analysing snowfall episodes began in Italy in 1681 with the publication of La Figura della neve (“The shape of snow”) by the Italian priest and natural philosopher Donato Rossetti (1633–1686). For a reliable knowledge of snow cover dynamics in the Samnite area, we generated a calibration dataset of 65 years in the timeframe 1683–2018 CE.

For the calibration purpose, a dataset of 97 years was compiled (88% in the period 1907–2018 and the remainder, discontinuously from 1683 to 1895, from historical data of the Benevento station), each year including the winter season starting in December of the previous year. They includes three groups of data, i.e., from the BNV, the MVOBS and the Abruzzi area (Table 1).

Table 1 Resources for model calibration and validation: sites/areas, geographical settings, data availability
Full size table

The BNV and the MVOBS host the longest DSG time-series in southern Italy. Snowfall data in the BNV (here seen as a point, not as an area) are organised by the Benevento city observatory and the Met European Research Observatory, 5-km distant from each other within the valley (Diodato 1997). The MVOBS, located at 1280 m a.s.l. and managed since 1884 by the Benedictine Community of Montevergine’s Abbey, is the oldest among the high-altitude meteorological observatories located in the southern Apennines (Capozzi and Budillon, 2017). The first 34 years are from the BNV, as derived from historical records, which are mostly mirroring extreme occurrences, i.e., years for which snow remained on the ground over many days. This is the case for the 16 years, which were characterised by snowfall events of remarkable impact for society and the agriculture (Diodato et al. 2019, 2020a, c, 2021). Snow cover duration was not obtained for other years, for which evidence of snowfall occurrence is not reported in written documents prior to 1973. More recently, the collection of 18 years of snow cover data in the BNV has been performed by the authors. In the second group of data, obtained from the MVOBS, the variation of snow cover records is larger, and the time-series of available data is more continuous because the station is snowy, and snowfall observations were performed with only occasional interruptions by the local religious community (SIMN 1921–1997). Here, the time-series of DSG is composed of 31 years. The third group of data is the result of a reconstruction for an area within the Abruzzi region in the 1939–1970 period, obtained by averaging snow cover data provided by SIMN (1921–1997) for three sites. Snow cover records at these sites have become less complete in more recent years, insufficient for reliable areal estimates. This study area also contains the Gran Sasso d’Italia Massif, with the Corno Grande peak (42° 28′ N, 13° 34′ E), the highest of the Apennines (2912 m a.s.l.), for which we do not have access to sufficient data of snow cover.

The above datasets were the basis for developing a semi-empirical model — Eq. (1) — relating DSG to weather and elevation inputs. Other datasets of observed DSG records or proxy data were derived for an independent validation of the model, within and beyond the CSA, covering a broad range of terrains (sites and areas) across the country (Table 1). Within the CSA, a qualitative validation was performed in the Samnite area, over the Matese Massif (~ 1800 m a.s.l.), by comparing model outputs with evidence of snow cover from a historical period (1781–1820) with the most recent phase (1979–2018). In order to evaluate the model beyond the CSA, a quantitative model assessment was performed over some complex terrains. They include Soveria Mannelli — locality San Tommaso (S. Tommaso OBS; SIMN 1921–1997) in southern Italy (south of ~ 39° N) and, in northern Italy (north of 44° N), the Parma observatory (Parma OBS; Diodato et al. 2020a), a northern Apennine area (Emilian Apennine, mean elevation ~ 1000 m a.s.l.; De Bellis et al. 2010), and the western alpine sector (mean elevation ~ 2000 m a.s.l.; Bovo et al. 2009).

At the Parma observatory (Parma OBS), centrally located in the Po River basin (44° 48′ N, 10° 19′ E, 49 m a.s.l.), snow cover observations are available since 1938 (Diodato et al. 2020a). Though occurring in the plain, snowfall at Parma reflects the important geographical and climatic complexity of this area, placed at the separation between the greatly snowy pre-Alps (to the north) and the Apennines (to the south), and the limited accumulations of snow to the east, where climate is less continental because of the proximity of the Adriatic Sea (Diodato et al. 2021). This makes the occurrence of snow events not easily predictable because mountain ranges create different orographic effects by hindering the movement of air masses (Enzi et al. 2014). In the same region (Emilia-Romagna), model estimates were also compared to the decreasing trend (from 1950) of the proxy-derived DSG records from De Bellis et al. (2010) for the northern Apennine (Emilian Apennine, mean elevation ~ 1000 m a.s.l.). Stretching around 44° 27′ N and 10° 47′ E, the Emilian Apennine is the northern sector of the Tuscan-Emilian Apennines on the border between the Emilia-Romagna and Tuscany regions. Likewise, we assessed model estimates across the western alpine sector (mean elevation ~ 2000 m a.s.l.; Bovo et al. 2009) and at the Soveria Mannelli — locality San Tommaso in the southern region of Calabria. In the Piedmont region, the mountainous area is formed as an arc that delimits its southern, western and northern borders, where the western Alps (or Piedmont Alps) are one of the most imposing mountain systems in Italy going from the Colle di Cadibona (44° 20′ N, 08° 21′ E, 436 m a.s.l.), a mountain pass delineating the boundary with the Apennine Mountains, to the Col Ferret (45° 53′ N, 07° 04′ E, 2537 m a.s.l.), an Alpine pass crossing the route to the Mont Blanc massif (with the highest point in western Europe at 4808 m a.s.l.). For the site of Soveria Mannelli (39° 05′ N, 16° 22′ E), located on hills below the Sila National Park, DSG data were obtained for 1936–1954, 1963–1972 and 1975–1985 from SIMN (1921–1997), which refer to the observatory of San Tommaso (S. Tommaso OBS, 820 m a.s.l.). The database of observed DSG in the BNV also contains the years 2019 and 2020 (with 3 and 1 days of snow cover, respectively), but these years were left out of the analysis because the CRU data are updated until 2018 (checked June 2020).

2.3 Development and design of the model

We adopted the criterion of parsimony following Mulligan and Wainwright (2013), stating that a parsimonious model is the one that accomplishes a desired level of explanation or predictive ability with the least number of parameters or processes. Based on this principle, we generated a relatively simple model that could be easily parameterised and validated, as well as reliably applicable to estimate DSG data over complex terrains in Italy. After Van der Knijff et al. (2000), northern and southern European hydrological conditions distinctly appear above 48° N and below 42° N, respectively, with less distinguishable conditions in the transition between these two limits. To increase predictability and minimize uncertainties in the modelling of annual DSG data, we built on a nonlinear multivariate regression with four input variables and six parameters, which were calibrated against the 97-year observational records available for the CAS dataset (the Snow cover data resources section). For this, we first acquired a comprehensive knowledge about potential drivers of DSG in both the Samnium and the Abruzzi, where snowfall occurs from December to April (more rarely in May), and then we used these factors (weather anomalies record, air temperature and its variability) to develop a parsimonious model for estimating DSG data. We could thus identify a sufficient number of yearly climatic explanatory variables to represent DSG dynamics, consistent with the sample of data available during the input selection process. First, we investigated the effects of single variables, or sets of variables, on DSG over some climatologically meaningful periods. Then, an iterative process (trial-and-error to compose relevant drivers) enabled us to explain the dynamics of DSG in relatively simple terms. A stepwise-selection logic (e.g., Talento and Terra 2013) was used to alternate between adding and removing terms. The following non-linear multivariate regression model was thus derived to estimate DSG (day year−1) over the CSA:

$$\mathrm{DSG}(\mathrm{CSA})=\mathrm{A}\bullet \left\{\alpha \bullet {{\left(1-AMV\right)}^{\upeta 1}\bullet \left(1+{SSI}_{l}\right)}^{\upeta 2}+\left[{\left(\beta -{T}_{w}\right)}^{\upeta 2}+({\gamma -{T}_{s})}^{\upeta 2}\right]\bullet \mathrm{VC}\right\}$$
(1)

where the AMV (°C) is the annual Atlantic multidecadal variability, SSI is the snow-severity index (l subscript standing for local), and Tw and Ts are the winter and spring mean temperatures, respectively (°C). The scale term α (°C−1), the position terms β (°C) and γ (°C), and the shape terms η1 and η2 are empirical parameters.

As DSG predictors, annually resolved time-series back to 800 CE are available for the two climate indicators AMV (Wang et al. 2017) and SSI (Diodato et al. 2019). The AMV is the fluctuation of North Atlantic Ocean surface temperatures on multi-decadal time-scales, expressing the Atlantic multidecadal oscillation (AMO) as part of its internal climate variability, and external forcing agents such as changes in solar radiation, volcanic eruptions and anthropogenic influences. It refers a warming/cooling pattern of North Atlantic sea-surface temperatures with a periodicity of ~ 60 years (e.g., Mazzarella and Scafetta 2012). The variability of North Atlantic Ocean temperatures exhibits an important influence on the climate of surrounding regions (e.g., O’Reilly et al. 2019), specifically acting on the Mediterranean hydroclimatic system in multiple ways (e.g., Maslova et al. 2017). Sea-surface temperature anomalies, as represented by the AMV, also induce atmospheric pressure gradients, redistributing air masses through the North Atlantic Oscillation (NAO) between subtropical and subpolar latitudes of the North Atlantic, and modulating the strength and latitudinal location of westerly flows (Hurrell 1995; Jones et al. 1997). The SSI is a grade of winter snow severity, developed for Italy (Diodato et al. 2019) on a range 0 (normal) to 4 (extraordinary), which contains a qualitative information related to snowfall anomalies at the local scale, as seen in historical documentation. Normal winter (SSI = 0) refers to long-term average snowfall or snowfall that went unnoticed because it was likely below the long-term average (this is also the case for years in which the chroniclers left no trace of snowfall in the records), without commenting on its severity or ramifications on society and the economy. At the other extreme, extraordinary snowy winter (SSI = 4) is characterised by very extreme sporadic events with a low century recurrence rate. Intermediate situations are represented by SSI = 1 (snowy winter), SSI = 2 (very snowy winter) and SSI = 3 (great snowy winter).

The limitations of the approach adopted to generate the Italian SSI time-series would be the excessive smoothing of local snowfall patterns. In order to capture the precipitation control on snow duration, which becomes dominant in high peaks, we based the extraction of local SSI inputs (SSIl) on the effective snowfall and written sources used by Diodato et al. (2019) to provide baseline index values for Italy (Supplementary material). The combination of AMV and SSIl is important to handle both temperature and precipitation controls on Italian snowfall. In fact, whether the AMV translates the impact of large-scale circulation anomalies on local temperatures, as altitudes increase it becomes progressively less likely for air temperatures to rise beyond the freezing point (Beniston et al. 2003).

Looking beyond the effects of large-scale temperature changes on snowfall, it appears that local seasonal temperatures are strongly related to the time during which snow remains on the land surface (Brown 2019). For the local temperature inputs, we used the seasonally resolved winter and spring temperatures (interpolated on a 0.5° by 0.5° grid), from Luterbacher et al. (2004) and updated from the climatic research unit (CRU) global climate dataset (Jones and Harris 2018). However, the regional reconstruction by Luterbacher et al. (2004) is known to underestimate the amplitude of low-frequency temperature changes, especially in winter and spring. In particular, it poorly reflects the variability of actual seasonal temperatures in southern Europe, where temperatures are far more variable than in central Europe (Diodato et al. 2013). Specifically, chronicles mark an unquestionable point of departure from the regional temperature reconstruction for the year 1684 in southern Italy (which also had the severest frost recorded in England; Manley 2011). For the exceptionally cold conditions of this year (with an expected return period of ~ 1000 years), we referred to the corrected mean seasonal temperatures derived by Diodato and Bellocchi (2011) for the BNV, i.e., 0.1 and 7 °C for winter and spring respectively (against 4.1 and 11 °C from the regional Luterbacher reconstruction). Adjusted temperature data were also derived for the north facing slope of the MM (at ~ 1800 m a.s.l., within a 0.5°-grid domain characterised by a sloping terrain at ~ 600 m a.s.l. on a average) by subtracting 7.2 and 8.4 °C (after Perosino and Zaccara 2006) to the regionally reconstructed values of Tw and Ts, respectively. We assumed vertical thermic gradients (lapse rates) of 0.48 °C 100 m−1 in winter and 0.56 °C 100 m−1 in spring (based on Rolland 2003).

A scale-transfer approach was required to extend the application of the model from single sites to areas. Since the spatial averaging is basically equivalent to the spatial integration (Chen et al. 1994), a viable possible scaling technique is averaging the site-scale model over an entire mountainside area by assuming that the object under study is a fractal (i.e., the site is a fractal element of the area). In this case, the scaling of the site can be reasonably modelled by a power-law process-area relationship, whose exponent η quantifies the spatial dependence (scaling nature) of the snow cover days (Zhang et al. 2013). This is an issue of determining how the climate of specific surface types or slopes differs from a grid-average climate, without formulating different model equations (Harvey 2013). In that way, a varying exponent not only provides a parsimonious description of the process under study, but it is also a generic mechanism governing the process. In particular, lower exponent values of mean areal response (compared to the response obtained at a single site) suggest that they are a construct (smoothing approximation) from individual sites (e.g., Spence and Mengistu 2019).

The term A is a scale function that balances the DSG with respect to a given elevation (E, m a.s.l.), modulated by position (ϕ), scale (σ) and shape (ψ) empirical parameters:

$$\mathrm{A}={\left[{\varphi }+\left(\sigma \cdot \frac{E}{{E}_{R}}\right)\right]}^{\uppsi }$$
(2)

The elevation ratio E/ER (after Diodato and Bellocchi 2020) reflects the proportion of DSG below or above a reference altitude (ER), which was set equal to 700 m a.s.l., i.e., roughly the average elevation of the CSA as well as the hilly and mountainous areas of Italy (337 m a.s.l. including plains, covering 23.2% of the Italian territory; CIRPIET 2004).

The term VC is a modified variation coefficient, that is, a standard deviation/mean ratio of winter–spring temperatures (T(w → s)), calculated across the current (t = 0) and previous (t =  − 1) year:

$$\mathrm{VC}=\left(\frac{\upeta \bullet \sqrt{SD}}{\mathrm{Mean}}\right)\left({T\left(w\to s\right)}_{t=-1}^{t=0}\right)$$
(3)

This is a key factor to infer seasonal temperature fluctuations as it tends to become larger with more persistent snow (Diodato et al. 2020a). This term modulates the temperature control across winter and spring, and acts as a memory effect of snow duration. We developed the hypothesis of Diodato and Bellocchi (2020) stating that site- or area-specific estimates of snowfall in peninsular Italy can be obtained from temperature variations in the BNV because this valley is located at the transition between central and southern Apennines and captures the temperature variations of the Italian peninsular sector. VC was thus calculated from BNV temperatures at all study-sites or areas of peninsular Italy (from the Po River Basin to the southern Apennines) while it was set equal to 1 in the alpine areas (i.e., at the highest latitudes and altitudes of Italy) where seasonal temperatures are mostly below 0 °C. Snowiness, not the variations of temperature, is the main cause of snow cover duration under these thermal conditions (Scaramellini and Bonardi 2001).

The concept of the model (Fig. 4) summarises the mechanisms mostly driving the common patterns of change of snow cover extent and duration (e.g., Brown et al. 2010; Brown 2019). It shows that in Italy, ~ 10–15% (central tendency) of the territory is covered by snow over wintertime, but in this season, the higher percentiles of snow cover extent deviate more from central values than in other seasons. This indicates that SSI > 0, which captures extreme snowy years (roughly above the median), is an important driver of snow cover as it reflects the winter snowfall regime. Snow cover responds to local temperature forcing over winter and spring months, while large-scale forcing (AMV) is expected to persist throughout the year. Compared to previous studies on snowfall modelling, limited to a few low-altitude sites (e.g., Diodato and Bellocchi 2020), the conceptual framework in Fig. 4 reflects the strategy adopted to extend the spatial generalisation capabilities of parsimonious modelling, which was necessary, given the geographical and climatic characteristics of the Italian territory (very dissimilar from North to South). This strategy, which includes elevation, and snowiness and teleconnection indices, allows a wide range of geographical and climatic variability to be captured and the model to be used to identify patterns of snowfall change in the CSA and elsewhere in Italy.

Fig. 4
figure 4

Monthly snow cover extent in Italy (50th percentiles in red, 80th percentiles in green, 98th percentiles in blue), with the associated driving factors. Percentiles (prcs) were calculated on data provided for the period 1966–2018 by Global Snow Lab – Rutgers University Climate Lab (https://climate.rutgers.edu/snowcover) via KNMI-Climate Explorer Climate Change Atlas (http://climexp.knmi.nl). The three main terms of Eq. (1) are reported

Full size image

The parameters of Eq. (1) were calibrated against the observed DSG dataset, optimized against the following criteria:

$$\left\{\begin{array}{c}\begin{array}{c}{R}^{2}=\mathrm{max}\\ \mathrm{MAE}=\mathrm{min}\end{array}\\ \left|b-1\right|=\mathrm{min}\end{array}\right.$$
(4)

R2 (ranging from 0, the worst, to 1, the optimum) is the goodness-of-fit of the regression between observed and estimated DSG data. MAE (from 0 day year−1, the best, to positive infinity, the worst) is the mean absolute error, the magnitude of which is a measure of the model error. The third condition approximates the unit slope (b = 1, optimum) of the straight line of the linear regression of the actual versus modelled data. We started with simple additive functions and then moved on to multiplicative functions of the four input variables (SSI, AMV, Tw, Ts). The physical interpretation of the studied phenomenon (Fig. 4) suggested the introduction of a multiplicative process. We iteratively added predictors, one-at-a-time, until modelling solutions with low MAE and high R2 were obtained. For the final selection, the third criterion, |b-1|, was considered. The calibration work was carried out through a trial-and-error (manual) process by comparing the model predictions with the observed data. The resulting solution was obtained through a series of > 50 forward and backward iterations and constitutes a satisfactory solution for all the three criteria (Diodato et al. 2020b). Then, we expanded the set of performance metrics used to evaluate the model. With the Pearson’s correlation coefficient r = √R2 (from − 1, the worst, to 1, the best), we measured the concomitant variation of estimates and observations. The normalized root mean squared error (NRMSE) is the root mean squared error divided by the mean of observed data, which facilitates the comparison between datasets with different scales. NRMSE values lower than 1 indicate satisfactory model accuracy. The Kling–Gupta index (− ∞ < KGE ≤ 1, Kling et al. 2012) was used as efficiency measure, with KGE > –0.41 indicating that the model is a better estimator than the means of observations (Knoben et al. 2019). We also calculated the Nash–Sutcliffe efficiency index (− ∞ < EF ≤ 1, optimal; Nash and Sutcliffe, 1970) as an indicator of model performance uncertainty, since values above 0.6 indicate limited model uncertainty, presumably associated with narrow parameter uncertainty (Lim et al. 2006). The analysis of variance (ANOVA) was applied to find out if all predictors were necessary (not redundant) for the modelling purpose.

Spreadsheet-based model development and assessment were performed with the analytical and graphical support of STATGRAPHICS (Nau 2005) and WESSA (Wessa 2009) statistical software. F-ratio P-values were used to present the statistical significance of the linear regression between estimated and actual data. The Durbin–Watson statistic (Durbin and Watson 1950, 1951) was used as a test for autocorrelation in the model residuals, considering that strong serial dependencies may induce spurious correlations (Granger et al. 2001). Exploratory and time-series analysis were carried out using AnClim (http://www.climahom.eu/software-solution/anclim).

3 Results and discussion

3.1 Model parameterization in the CSA

The calibrated parameter values of Eq. (1) are α = 3.45, β = 9.0 °C, γ = 15.0 °C, η1 = 2.0 and η2 = 2 for single stations and 1.7 for the surface area of a terrain. For Eq. (2), calibrated parameter values are: φ = 0.35, σ = 0.50 and ψ = 0.83. There is statistically significant relationship between observed and modelled values (P < 0.05). The value of the R2 statistic of observed versus modelled data (Fig. 5a) indicates that the fitted model explains 88% of the observed variability. Only a few data points (6%) are outside the bound of 95% prediction limits for new observations (light pink in Fig. 5a). The Nash–Sutcliffe efficiency value (EF = 89) indicates limited uncertainty in the model estimates. The mean absolute error (MAE), used to quantify the amount of error, was equal to 6.0 day year−1, which is lower than the standard error of the estimates (8.8 day year−1). Other performance metrics (NMRSE near 0 and KGE near 1) also reveal that model estimates are satisfactory. The Durbin–Watson statistic 1.91 (P = 0.30) indicates that there is no significant serial autocorrelation in the residuals.

Fig. 5
figure 5

Scatterplot of observed (y) and modelled (x) DSG for the calibration stage, with a regression (y = 0.73 + 0.99·x; black line) and 1:1 line (y = x; bold pink line) with the inner bounds showing 90% confidence limits (dark pink area), and the outer bounds showing 95% prediction limits for new observations (light pink area), b histogram of residuals, and c Q-Q plot with sample versus theoretical quantile values. Performance metrics are reported in a

Full size image

The pattern of residuals approaches a normal distribution (Fig. 5b; normality test, P > 0.05, Jarque and Bera, 1981). The Q-Q plot (Fig. 5c) exhibits a distribution of sample-quantiles around the theoretical line (only slightly tailed).

Taking the four input variables individually, we note that winter and spring temperatures are negatively correlated with DSG (r =  − 0.65 and r =  − 0.64, respectively), which is little or weakly correlated with AMV (r =  − 0.20) and SSI (r = 0.48). The four input variables are mostly uncorrelated: AMV-SSI r =  − 0.14, AMV-Tw r =  − 0.01, AMV-Ts r = 0.05, SSI-Tw r = 0.23 and SSI-Ts r ~ 0.00. To determine whether \(\mathrm{DSG}(\mathrm{CSA})\) can be simplified, we fitted a multiple linear regression between DSG and two independent, uncorrelated (r = 0.01) terms: \(\alpha \bullet {{\left(1-AMV\right)}^{\upeta 1}\bullet \left(1+{SSI}_{l}\right)}^{\upeta 2}\) and \(\left[{\left(\beta -{T}_{w}\right)}^{\upeta 2}+({\gamma -{T}_{s})}^{\upeta 2}\right]\bullet \mathrm{VC}\) of Eq. (1). Both variance components of the model are highly significant (P ~ 0.00 for both terms) and the model cannot be simplified further, resulting in a stable, interpretable and generally useful model (Royston and Sauerbrei 2008).

3.2 Qualitative validation (within the CSA)

We provide a qualitative assessment of model accuracy within the CSA, in the absence of independent sets of snow cover observations. This qualitative validation was performed by comparing the modelled evolution of snow cover in the Matese Massif (at ~ 1800 m a.s.l.) in the past four decades (1979–2018; Fig. 6b) with earlier estimates in a period of four decades (1781–1820; Fig. 6a) of the late LIA, which roughly overlaps with the Dalton solar minimum (~ 1790–1830 CE), a period of below-average global temperatures characterised by reduced solar activity and enhanced volcanic events (Wagner and Zorita 2005). Four decades give us a sufficient time span to compare patterns of precipitation variability (Perini 1990). An accumulation of permanent snow was typical in the highest mountains surrounding the BNV during the late LIA, when summer snowfields used to whiten the peaks of these mountains (including the MM), as referred by Quattromani (1827, p. 18): “In places of average latitude the cold is more intense for the particular localities to which the Apennines contribute, part of which is covered with snow in the winter, and in the same summer the snows do not clear the tops of Corno, Matese and the Majella Mountains.”

Fig. 6
figure 6

DSG time-series (blue curve) in the Matese Massif (~ 1800 a.s.l.) as estimated with Eq. (1) for the periods a 1781–1820 and b 1979–2018, and c monthly (from January to May) snow-depth observations (cm) in 2012 at Campitello Matese ski resort (41° 27’ N, 14° 23’ E, 1400 m a.s.l.), as from https://www.skiinfo.it/molise/campitello-matese/storico-neve.html?y=0&q=top. In graphs a and b, pink bands represent the summer period (from June 1 to August 31)

Full size image

Figure 6a (blue curve) shows that the north face of MM peaks (over ~ 1600–2000 m a.s.l.) remained snowy during the summer months of June, July and August (Fig. 6a, pink horizontal band) over 70% of the years 1781–1820 (a period consistent with the historical description provided by Quattromani 1827). This is also supported by the testimony offered by Tenore (1827, p. 44), who wrote “In the Matese valleys the snow remains most of the year”, thus suggesting that the highest MM peaks were so cold in the early nineteenth century that snow did not melt in summer. At that time, snowfield was prominent in the MM: “Our king of the mountains the high Matese, to whom the snows freeze, even when in Leo the Sun dwells, covers the chin, and the white head” (Vv.Aa., 1873, p. 510). It had to dazzle the travellers who walked along the Seapinum plateau, located at ~ 500 m a.s.l.: The wide and vast dewlap of Matese, then, was truly majestic: all the mountains were in sight, and above them emerged superb, drawing in the blue background of the air, the beautiful profile of its round dome, Mount Miletto (Vv.Aa., 1892, p. 226).

While pre-LIA conditions indicate (e.g., Ivy-Ochs et al. 2009) smaller glaciers than today during Roman times (until ~ 200 CE) and during the central part of the Middle Ages (roughly from the ninth century to the thirteenth century), the glacial tongues of today can be considered as the last legacy of the past uniform snowfield of great extent during the LIA that in some years survive at the beginning of summer (https://www.meteoinmolise.com/nevaio-in-molise-le-straordinarie-foto-di-martusciello-fino-a-20-metri/34409). A comparison of the perspective view of the north face of Mount Miletto (as shown in Fig. 7a) with a recent picture taken from Campitello Matese ski resort (Fig. 7b) shows the areas that were occupied by the historical (permanent) snowfield (a LIA phenomenon of when Italian glaciers likely reached their Holocene maximum advance, e.g., Orombelli 2011), which has been shrinking over time. They appear as shadowed silhouettes on the concave side of the mountain (Fig. 7b). This glacier history of Italy broadly follows the one for the Alps (Holzhauser et al. 2005) and for Europe in general (Solomina et al. 2016).

Fig. 7
figure 7

a Perspective view of the north side of Matese Massif (with Mount Miletto, 2050 m a.s.l.) from the Saepinum Plateau (548 m a.s.l.), b picture of the massif from Campitello Matese (41° 27′ N, 14° 23′ E, 1400 m a.s.l.) ski resort (http://www.campitellomatese.org/it/dove-siamo/il-comprensorio-di-campitello-matese-009.html). Red rectangle in a locates the historical snowfield highlighted in b

Full size image

This contrasts with the most recent warm decades. Figure 6b shows that summer snow cover in the MM peaks is absent in most of the timeframe 1979–2018 (pink horizontal band), with the exception of a few years (~ 25%) for which snow cover was modelled in early summer. With 195 days of snow cover (estimated from 1 December 1992), we identified 1993 as the last year in which snow covered the MM peaks until the beginning of June. The year 2012 was particularly snowy in that area (Bisci et al. 2012), but snow on the ground likely disappeared in May (Fig. 6c, red square). Our estimates are in agreement with snow-depth monthly records collected in 2012 at Campitello Matese (a ski resort of the MM, located at 1450 m a.s.l.), indicating that snow only stayed on the ground until the beginning of May.

The progressive withdrawal of the snow cover, as simulated by our model in the MM (Fig. 6b), is also in agreement with the glacier contraction that occurred in the adjacent Abruzzi Apennines. In fact, only traces remain today (Pecci et al. 2008) of the Calderone glacier (42° 28′ N, 13° 33′ E, ~ 2900 m a.s.l., Fig. 8b), located in the Gran Sasso d’Italia massif (Fig. 8a). With our qualitative validation, we claim that the model can be used to represent interdecadal and long-term trends without speculating on its use to estimate single-year DSG data. In the next section, this qualitative validation (performed within the same area where the model was calibrated) is complemented with a quantitative validation to ensure that the model complies with its intended purpose, i.e., to represent snow cover dynamics in Italy beyond the calibration area.

Fig. 8
figure 8

a Landscape of the Gran Sasso D’Italia mountain group (https://www.touringclub.it/itinerari-e-weekend/cinque-idee-per-andare-in-montagna-nel-ponte-del-25-aprile/immagine/9), and b north face of Corno Grande with traces of the southernmost Calderone glacier in July 2019 (https://www.caputfrigoris.it/calderone.htm)

Full size image

3.3 Quantitative validation (beyond the CSA)

Figure 9 compares DSG predictions to observed trends at Parma OBS (a), Emilian Apennine (b), Piedmont Alps (c) and S. Tommaso OBS (d). Performance metrics (r > 0.6, NMRSE < 0.5 and KGE > 0.5) are satisfactory at each site or area. The trend towards anomalously warm conditions in the last decades indicates that increasing temperatures represent the main factor triggering the decline of snow-duration days on the ground, especially in northern Italy. The satisfactory performance obtained at Parma OBS (EF = 0.81) is remarkable because it overcomes two major challenges for DSG modelling.

Fig. 9
figure 9

Comparison of modelled (orange curve) and observed (black curve) DSG over time at validation sites and areas (Fig. 1c): a Parma OBS (1938–2002), b Emilian Apennine (1950–2002), c Piedmont Alps (1932–2002), and d S. Tommaso OBS (1936–1954, 1963–1972, 1975–1985). Elevation above sea level (ele., m) is reported for each site, while the average elevation (avg. ele., m) is given for areas. Performance metrics of model goodness are described in section Development and design of the model

Full size image

First, this indicates that the model can be employed to investigate changes in the snowfall dynamics of northern Italy’s lowlands. Here, snowfall episodes are mostly caused by the softening of the cold air layer (“cushion”) forming at low altitudes. Under such temperature inversion conditions, snow begins to fall in the plains as the moist air (warmer and lighter) that arrives with perturbations causes the gradual warming of the underlying layers until dismantling the pockets of cold air that are there. This phenomenon is more common in the western sector of the Po River Basin (Steiner et al. 2003), which has a distinctly continental climate and where the distribution of the reliefs favours the stagnation of cold air in the valley (Caserini et al. 2017). To simulate DSG in lowlands, the model can be informed of this higher layer’s control on the lower layer by setting as input a higher elevation level than the plain level. The latter is 49 m a.s.l. for the Parma OBS, for which we obtained the best DSG estimates setting E = 180 m a.s.l. This elevation can thus be used as the minimum elevation threshold for the application of the model, though its value could be revised in the future to accommodate other datasets.

Second, the model offers robustness against anthropogenic disturbance in urban environments. The latter is relevant for Parma, where Zanella (1976) showed an average difference of 1.4 °C between urban and rural temperature data in the period 1959–1973, and reported that the difference varied seasonally (especially in spring and summer). Heat islands may thus have been a contributing factor of the snow cover decline (e.g., Musco 2016) observed in this urban site (from about 80 day year−1 in the 1940s to < 20 day year−1 in the most recent years), which appears to be sufficiently captured by the modelled time-series.

An important decline of DSG values was also observed in the northern Apennines (Fig. 9b; EF = 0.61), where 100 day year−1 with snow on the ground have no longer been exceeded over the past decades, starting in the 1960s. This decline was interpreted by the model, which can thus be used to support studies on mountain ecosystems, which are sensitive and reliable indicators of climate change. In fact, Petriccione and Bricca (2019) report that snow cover contraction and an increase in mean and minimum annual temperatures could have resulted in as much as a 50–80% change in plant-species composition with respect to the total number of species observed since the 1980s in the central Apennines. In turn, this could lead, in the medium- or long-term, to degradation processes affecting high-elevation plant communities of the Apennines (including the local extinction of most of the cold-adapted species), due to their low resilience. While the performance of the model is satisfactory in southern Italy as well (Fig. 9d; EF = 0.86), lower correlation and efficiency (EF = 0.38) values obtained in the Piedmont Alps (Fig. 9c) may indicate that such high altitudes mark the upper limit of applicability of the model (where it can be implemented in support of geomorphological process studies and hazard-risk mapping, e.g., Keiler et al. 2010).

4 Conclusions

Recent research has highlighted the difficulty of monitoring the length of snow on the ground in complex terrains. Although satellite measurement data have been used to detect the snow cover duration in broad mountain systems, numerous questions and uncertainties remain to be explored. The three main challenges in resolving these uncertainties are (1) improving the knowledge of the main factors driving snow cover dynamics; (2) promoting the understanding of the strengths and limitations of modelling approaches when making climate interpretations of such dynamics and (3) better incorporating and accounting for different scales of observation and uncertainty in the analyses of changing snow cover durations. New satellite missions, and better access to heritage data (documentary data), are improving our ability to monitor remotely sensed snow cover changes. However, the generation of new such data alone will not provide solutions to many of the long-standing conceptual challenges. For that, we need to look forward by quantifying contemporary and future snow cover duration changes at the same time as we look back through the analysis of historical snow cover data. In general, the use of data deriving from a variety of sources is not straightforward in reconstructing snow cover durations. The main novelty of this study is the development (with the support of publically available resources) of a relatively simple model to reconstruct the variability of DSG at the country level, based on historical and recent data. With the modelling of annually resolved DSG time-series from Italy, we have deepened our understanding of the dynamics of snow cover in complex terrains and reported the consistency of our estimates with instrumental data and historical evidence. We report the results of our study with respect to the main drivers of the observed changes in DSG data. They include large-scale changes in the modes of climate (i.e., the Atlantic multidecadal variability), identifying in this way the need for future research on the subject (e.g., on the non-linear dependence of snow cover duration from temperature). Hence, the datasets gathered and the model developed in this study offer us the possibility to make statements about the relationship between DSG in Italy and the wider regional climate which, to our knowledge, are unprecedented in their detail and confidence. However, though our study has revealed that DSG patterns may be represented by a combination of local and large-scale climate factors, their relationship with snow cover duration mechanisms in Italy should be regarded as still tentative, and in need of confirmation by more extensive studies. In the meantime, we hope that this study will contribute to enhanced knowledge of the snow cover dynamics in the Mediterranean areas where they are still poorly known, and where the approach could be suitable for applications, provided that the model parameters are documented for territories other than Italy.