Multi-Decadal Variability in the Snow-Cover Reconstruction at Parma Observatory (Northern Italy, 1681–2018 CE)

Diodato, Nazzareno; Bertolin, Chiara; Bellocchi, Gianni

doi:10.3389/feart.2020.561148

ORIGINAL RESEARCH article

Front. Earth Sci., 30 October 2020

Sec. Hydrosphere

Volume 8 - 2020 | https://doi.org/10.3389/feart.2020.561148

This article is part of the Research Topic Advanced Time Series Analysis in Geosciences View all 13 articles

Multi-Decadal Variability in the Snow-Cover Reconstruction at Parma Observatory (Northern Italy, 1681–2018 CE)

Nazzareno Diodato¹

Chiara Bertolin²^*

Gianni Bellocchi^1,3

¹Met European Research Observatory – International Affiliates Program of the University Corporation for Atmospheric Research, Benevento, Italy
²Department of Mechanical and Industrial Engineering, Faculty of Engineering, Norwegian University of Science and Technology, Trondheim, Norway
³Université Clermont Auvergne, INRAE, Clermont-Ferrand, France

Emerging negative trends in snow depth and cover days highlight the challenges posed by changing snow patterns around the world. They suggest that snow-dependent regions in southern Europe could be affected by these changes because the number of days with snow on the ground (DSG) determines soil processes and water-flow in rivers, streams, lakes and reservoirs. We present here the first homogeneous, annually-resolved (from October to April), multi-centennial (1681–2018 CE) DSG time-series for the Parma meteorological observatory (OBS), in northern Italy, which to date is also the longest DSG series reconstructed in the world. DSG data are in fact still poorly documented and misunderstood due to the limited and fragmentary data measurements of the past. DSG recording only began in 1938 at Parma OBS. To generate the long-term annual DSG time-series at the study site, we develop a model consistent with calibration (1938–1990 CE) and validation (1991–2018 CE) samples of observed data. We show that the variability of DSG depends on winter precipitation and air temperature, as well as on winter-spring temperature variability, suggesting that long sequences of DSG are dominated by cold air masses in years with cold weather and high variability. Modeled DSG data show a downward trend from the 19th century, in the transition period from the cold of the Little Ice Age to the warmth of modern times, followed by a more rapid decline in the five most recent decades. The DSG at Parma OBS appear to have followed over the last century trends similar to those observed throughout Eurasia and across the Northern Hemisphere, where a marked decline of snow-cover duration has been reported in the transition seasons (spring and autumn).

Introduction

Publio Virgilio Marone – Georgiche, I century B.C., Libro III (Italian translation from Latin: Clemente Bondi, 1801; our translation to English).

Every snowfall deposits its own layer on the soil surface. The timing, duration and persistence of the snowfall, which depend on climatic conditions, influence the amount of water stored in soils and its use for agriculture and other sectors such as industry, urban and rural domestic water supply (e.g., Lute et al., 2015; Barnhart et al., 2016; Wang et al., 2018a; Jennings et al., 2018). Large-scale snow-cover anomalies also cause important changes in the diabatic heating of the Earth’s surface by enhancing the fraction of solar radiation reflected away by the surface (i.e., the surface albedo), thus becoming an essential component of the terrestrial radiation balance (e.g., Cohen and Rind, 1991; Goward, 2005; Sandells and Flocco, 2014). Snow cover affects the timing and magnitude of flow peaks generated by snow melting (Wang et al., 2006). By delaying the transfer of precipitation to surface runoff and infiltration in catchment areas, it also influences flow reduction and the extent of the flow network during summer baseflow (Yarnell et al., 2010; Godsey et al., 2014).

Knowledge of past climate is an important key to understanding long-term snow-cover variability. For instance, capturing anomalously cold and snowy winters can be consistent with, and help explain, a persistence of snow cover on the surface in late spring (after Jungclaus, 2009; Handmer et al., 2012; Enzi et al., 2014). In particular, the interannual variability of the number of days of snow on the ground (DSG) is an important part of the climate signal to detect potential changes of snowfall intensity and spatial distribution of snowpack variables that may have important impacts on both the environment and the society (Changnon and Changnon, 2006; O’Gorman, 2014). Our understanding of the DSG characteristics in several regions of the world is limited because of the short records available and a poor knowledge of the complex weather and climate patterns that occur locally. In Europe, the snowfall dynamics depend on the latitude and atmospheric circulation patterns, modulated by local orographic situations (Croce et al., 2018; Kretschmer et al., 2018). Polar maritime air masses originating from over the Atlantic collide here with the polar continental air masses connected with the Asian high pressure (Bednorz, 2004). The temperature control is dominant in the transient snow regions where the mean winter temperature is slightly below, and often crosses, the melting point (Diodato and Bellocchi, 2020). The situation is different in regions where the mean winter temperature is well below 0°C, as is the case of the Alps, where increases in DSG are controlled mainly by precipitation inputs (Clark et al., 1999). At lower-elevations, DSG are instead affected principally by air temperature and only secondarily by precipitation conditions.

Data on snowfall and the persistence of snow cover on the surface are becoming increasingly important due to the high variability of snowfall rates worldwide (e.g., Davi et al., 2012), and the impact that snow cover can have on society, agriculture and water resources. A quantitative assessment of the long-term and interannual DSG variability is noticeable for the identification and framing of signals of climate change, the validation of climate models, and a better understanding of the interactions among the different spheres of the Earth system, including the geosphere, biosphere, atmosphere, hydrosphere and cryosphere (Brown, 2014). Documentary sources, usually in manuscripts or annotations in different formats, provide evidence and useful information to study the variability of the climate over historical periods (Glaser and Riemann, 2009; Dobrovolný et al., 2010). However, the quality and quantity of sources are often unevenly distributed in space and time (Bradley and Jones, 1992; Mann et al., 2000; Brázdil et al., 2005; Adamson, 2015). This is especially true for snowfall, considering that in situ continuous snow monitoring is an arduous task (De Walle and Rango, 2008) and satellite remote sensing data capturing the snow-cover evolution are limited to recent decades (e.g., Pimentel et al., 2017). Also for historical times, snow characteristics data remain still poorly documented and understood due to limited and fragmented measurements (Kunkel et al., 2016; Figure 1). One of the longest-running snowfall records is the series of Parma observatory (Parma OBS), in the Padan plain (northern Italy), where however recording of DSG only started in the winter season 1937–1938. This coincides with the setup of new measuring meteorological facilities at ground level and no longer on the tower of the OBS, the suburbs and the city center.

FIGURE 1

FIGURE 1. (A) The Parma meteorological observatory from a historical picture; (B) monthly summary of precipitation records (rainfall and snowfall) for the year 1839 (note the absence of recording in days with snow on the ground). Historical pictures are: (A) made freely available from Fabrizio Bonoli (University of Bologna, Italy) through the web-based content store (in Italian) of Ravenna Planetarium (http://planet.racine.ra.it/testi/astr_it.htm) [Accessed September 10, 2020]; (B) an excerpt of a printed image from Colla. (1840).

Centrally located in the Po River Basin (PRB), the Parma OBS started weather observations on 1777 thanks to the local Jesuit community with the support of the University of Parma, but only in the 20th century regular and continuous daily snowfall measurements were recorded. Diodato et al. (2020a) provided the series of monthly snowfall reconstruction in Parma over the 1777–2018 period, with the detailed analysis of metadata. In this regard, Italian regions hold among the world’s longest monthly snowfall time series (Enzi et al., 2014), going back to 1780 at Rome (Mangianti and Beltrano, 1993), 1788 at Turin (Leporati and Mercalli, 1994) and 1884 at Montevergine (Diodato, 1997). Nevertheless, Parma OBS possesses the longest continuous record, as it goes back to 1777, although we know (Camuffo and Bertolin, 2012) that meteorological observers were already active in Parma from 1654 to 1660, supported by the Grand Duke of Tuscany Ferdinand II de’ Medici (1610–1670).

The current work presents the first annually-resolved reconstruction of DSG since 1681 (i.e., the snow winter season going from October 1680 to April 1681) for the Parma OBS. A day with snow cover is a day with snow of at least 0.01 m depth observed in over 50% of an open neighboring 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 for the purpose of this study (e.g., Falarz, 2004; Petkova et al., 2004). We first acquired a comprehensive knowledge about potential drivers of DSG in the PRB (Materials and Methods), where snowfall occurs from October to April, and then we used these factors (amounts of precipitation and air temperature) to develop a tight or “parsimonious” model for reconstructing annual DSG data over the period 1681–2018 CE (Results and Discussion). This allowed us to capture a wide range of climate variability and identify patterns of snowfall changes reported and discussed in section.

Materials and Methods

Environmental Setting and Data

Parma OBS (44° 48′ N; 10° 19′ E, 49 m a.s.l.) is located in the central part of the PRB, in the Italian administrative region Emilia-Romagna (Figure 2A). This area is surrounded by the pre-Alps to the north and by the Emiliano Apennine to the south (Figure 2B). The interaction between these morphology and weather characteristics governs the occurrence of DSG (whose measurements are available since July 1937) under different conditions. The Parma OBS (black dot in Figure 2C) outlines a separation between the snowy peaks over the northern and western pre-alpine and alpine chains and the Apennines to the south where the climate is less continental. Complex reliefs, different slope exposures, as well as various distances to sea (Figure 2B), exert highly contrasting effects on snow-cover depth, persistence and spatial extent. The colored bands in Figure 2C do explain the amplitude of sub-regional variations around the Parma OBS. For instance, south of Parma the average snow-cover duration rises up to ∼80 days year⁻¹ in the Apennines.

FIGURE 2

FIGURE 2. Geographical setting: (A) North Italy (blue square); (B) the Po River Basin (black curve) with main locations and Parma OBS; (C) related winter (November-March, 2000–2011) snow-cover duration. Maps are authors’ own elaboration from free, public domain images: (A) ESRI (http://www.esri.com) [Accessed September 10, 2020]; (B) Elastic Terrain Map (http://elasticterrain.xyz) [Accessed September 10, 2020]; (C)Dietz et al. (2012).

During the year, cyclonic Atlantic air masses and areas of anti-cyclonic airflows from central and eastern Europe alternate. Fields of variable winds are associated with wet air masses from southwest, and dry and continental winds from northeast. Snowfall increases with elevation until 2000 m a.s.l. Figure 3 (gray bars) refers to observed monthly snow cover days. At Parma OBS, which can be assumed sufficiently representative of the PRB area, the average of DSG in the period 1938–2018 is about 20 days year⁻¹. Over this period of 81 years, the month of January had the greatest number of days with snow cover, with the soil covered by snow for about 10 days on average. On February, the number of days with snow on the ground roughly doubled the days of snow coverage in December. March and November had a similar number of days with snow cover, while the lowest number was registered in October. A similar distribution is noticeable for the 98th percentiles (Figure 3, empty bars).

FIGURE 3

FIGURE 3. Monthly mean distribution of days of snow on the ground (DSG) at Parma OBS (gray bars) with the related 98th percentile (empty bars), calculated over the period 1938–2018.

The most abundant snowfall in the Alps and Apennines occurs during cold periods, in particular with the arrival of warmer and moist air from the south, with a rise in temperatures up to about 0°C. If plain areas are exposed to a degree of cold sufficiently intense, with temperature values below 0°C, snowfall and its persistence on the ground are also possible at very low altitudes (Bettini, 2016). This is what happens in the central-western sector of the PRB, when warm and humid air coming from the Mediterranean Sea slips on the cushion of cold air previously trapped in the low layers of mountain landscapes in the Alps and Apennines. There is a chance that the arrival of below-freezing temperatures can also bring snow along the PRB. In the presence of air masses close to saturation, very cold air suddenly entering the PRB raises near saturation leading wet air to condensation. This phenomenon, although rare, can be observed in the PRB when cold Siberian air penetrates westward into northern Italy, while very damp air or even a misty blanket stagnates in the plains. In these cases, snowfall is weak, of short duration and limited to the plains or the lowest slopes of the pre-Alps.

As covariates for the modeling of DSG, we referred to winter precipitation, and winter and spring air temperatures. Winter and spring are defined as periods from December to February, and from March to May, respectively. For the precipitation input, we used the seasonally-resolved data with 0.5° resolution, as arranged from Pauling et al. (2003). Likewise, we used seasonally-resolved winter and spring temperatures as arranged by Luterbacher et al. (2004) since 1500, and interpolated at the Parma OBS grid point. Both datasets were updated from the Climatic Research Unit global climate dataset (Jones and Harris, 2018). The database of observed DSG contains also the year 2019 (with 4 days of snow cover on the ground) but this year was left out of the analysis because the Climatic Research Unit data are updated until 2018 (checked on June 2020). For both winter precipitation and temperature, the analysis of the yearly data used (1681–2018 CE) shows that they are compliant with a normal distribution using the Kolmogorov-Smirnov test p > 0.05 (Supplementary Appendix A).

Development and Parameterization of the Statistical Model

The principle of parsimony as described in Mulligan and Wainwright (2013) states that a parsimonious model is the one with the greatest explanation or predictive power and least parameters or process complexity. Based on this principle, we developed a parsimonious model that could be easy parameterized and validated, and is reliable and applicable to the reconstruction of long-term historical DSG data. To increase predictability and minimize uncertainties in the modeling of annual DSG data, we built a nonlinear-multivariate regression with three input variables and three parameters considering the original observational 1938–2018 CE timeframe. The data resource of annual DSG data (81 years) was separated into two sub-sets to use roughly two-thirds (53 years) for calibration (1938–1990 CE) and one-third (28 years) for validation (1991–2018 CE). The calibration effort was thus anchored to the reality of snow-cover duration in the past to enable a reasonably accurate reconstruction of historical snow-cover days. A number of monthly climatic explanatory variables was considered during the input selection process. First, in order to reduce the number of inputs, 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 decompose relevant drivers) enabled us to explain the dynamics of DSG in relatively simple terms. A stepwise selection logic was used to alternate between adding and removing terms. The following nonlinear multivariate regression model - DSG(Parma) - was thus derived to estimate DSG (day year⁻¹) at Parma OBS:

D S G (P a r m a) = A \cdot [l n (P_{W}) \cdot {(β - T_{w})}^{η} \cdot V C (\frac{S D^{η}}{M e a n}) (T {(w \to s)}_{t = - 1}^{t = 0})] (1)

where: A (day year⁻¹) is a scale parameter corresponding to the number of DSG when the term between brackets is equal to unity; β (°C) and η are process parameters. Winter temperature, T_w (°C), is negatively related to DSG, as mediated by the shape in parameter η. Winter (w) and spring (s) temperature variability accounted for by the term VC is indicative of the instability of air masses near the Parma OBS. This approach discloses the association between the snowfall response variable DSG (Parma) and multiple predictors, such as the winter precipitation amount (P_w), winter mean air temperature (T_w), and winter-spring modified variation coefficient VC, that is, a standard deviation/mean ratio of winter and spring temperatures (T(w→s)), calculated across the current (t = 0) and previous (t = −1) years. This variability term is a key factor to infer inter-seasonal temperature fluctuations, as it is larger in the presence of snow.

The concept of the model (Figure 4) summarizes 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 ∼10–15% of the Parma territory (central tendency) is covered by snow over wintertime. In this context, cool-season precipitation amounts are particularly important because the majority of annual precipitation occurs during the wintertime at medium-high latitudes, and preferentially 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). The logarithmic transformation of the total (solid and liquid) winter precipitation (ln(P_w)), which acts as snow-depth index in the model, was applied to underweight the amount of freshly fallen snow associated with broad-scale mixed precipitation over flat areas (e.g., Meir et al., 2016; Fehlmann et al., 2018). Snow-cover response to temperature includes the spatial extent (0.5°-grid resolution of winter temperatures) and the temporal scale (winter and spring months) over which the thermal forcing becomes important, 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 modeled by a power-law process 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 a specific surface type differs from a grid-average climate, without formulating different model equations (Harvey, 2013). In particular, the η exponent does not only provide a parsimonious description of the process under study, but it is also a generic mechanism governing the process. In particular, higher exponent values obtained at a single site (compared to mean areal response) suggest that they can be regarded as a method to downscale areal (smoothed out) approximations to finer units (e.g., Spence and Mengistu, 2019).

FIGURE 4

FIGURE 4. Monthly snow-cover extent in Parma territory (100-km square grid cell), with the associated driving factors. Monthly mean values (red line) 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) [Accessed September 10, 2020]. The three main terms of Eq. 1 are reported.

We have also found that a strong correlation (r = 0.84) exists between snowfall frequency and snow-cover duration at Parma in the period 1938–2018 CE (Supplementary Appendix B). However, since snowfall frequency measurements (i.e., snow days per year) only started in 1777 at Parma OBS, we excluded this input from the model, which was built on the temperature- and precipitation-based proxies allowing to extend the reconstruction period back to 1681. With this long-term series, we can detect patterns of climate change explaining fluctuations and tendencies in DSG data. Going back about 100 years before the beginning of the snowfall-frequency observations at Parma OBS, we could in fact start our DSG series at the turning point of the 17th and 18th centuries, at the so-called Maunder Minimum, which started in 1645 and lasted until 1715, a period characterized by exceedingly rare sunspots and lower-than-average temperatures (Eddy, 1976).

Spreadsheet-based model development and assessment were performed with the analytical and graphical support of STATGRAPHICS (Nau, 2005) and WESSA (Wessa, 2009) statistical software. The Mean Absolute Error (MAE, day year⁻¹) was calculated to quantify the differences between actual and modeled DSG values, while the Kling-Gupta index (−∞<KGE ≤ 1; Kling et al., 2012) was used as efficiency measure (with KGE > –0.41 indicating that a model is better performing that the mean of observations as a benchmark predictor; Knoben et al., 2019). The Nash-Sutcliffe efficiency (−∞<EF ≤ 1, optimum; Nash and Sutcliffe, 1970) was also calculated as an uncertainty indicator of the model performance because greater values than 0.6 indicate limited model uncertainty, likely associated with narrow parameter uncertainty (Lim et al., 2006). With the determination coefficient (0 ≤ R² ≤ 1, optimum), or the correlation coefficient ( $r = \sqrt{R^{2}}$ ), and the slope of the regression actual vs. modeled data (b = 1, optimum), we selected the set of important covariates for the parsimonious model for estimating DSG (after Mulligan and Wainwright, 2013). F-ratio p-values were used to present the statistical significance of the linear regression between actual and estimated data, and of the inputs’ relationship to the dependent variable. The Durbin-Watson statistic (Durbin and Watson, 1950; Durbin and Watson, 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). A MATLAB toolbox (https://noc.ac.uk/business/marine-data-products/cross-wavelet-wavelet-coherence-toolbox-matlab) was used for wavelet coherence analysis (Grinsted et al., 2004).

Results and Discussion

Model Parameterization and Evaluation

For the reconstruction of a long-term annually-resolved DSG series, we first calibrated the DSG(Parma) model. The calibration work was performed through a trial-and-error process comparing the model estimates with observations until small MAE and large R² values were obtained. Then, for the final selection of the parameter values, a third criterion (KGE closer to 1) was additionally involved. The calibrated parameters are: A = 0.0404 day year⁻¹, β = 6.80°C, and η = 2.00. The value of the R² statistic of observed (y) vs. modeled (x) data (Figure 5A) indicates that the fitted model explains 72% of the observed variability. The regression line y = 4.36 (±8.39) + 0.83 (±0.27)·x has intercept (a = 4.36) and slope (b = 0.83) with a relatively high standard error of the intercept (±8.39 days year⁻¹) and a root mean standard error of the fit equals to 0.24 day year⁻¹. This indicates the model’s lesser predictive ability for near-zero DSG, there is a superior predictive ability of Eq. 1 compared with a purely linear multivariate regression against the same independent variables, which produces several (unrealistic) negative values of DSG in both the calibration and validation periods (Supplementary Appendix C). The Nash-Sutcliffe efficiency value obtained in the calibration stage (EF = 0.73) indicates limited uncertainty in model estimates. Since the F-ratio p-value is less than 0.05, the linear regression between actual and estimated data is statistically significant. The mean absolute error (MAE), used to quantify the amount of error, was equal to 8.9 days year⁻¹, which is lower than the standard error of the estimates (10.9 days year⁻¹), while the Kling-Gupta index of 0.6 suggests a model of sufficient quality. The Durbin-Watson statistic (DW = 1.72, p = 0.14) reveals that there is no significant autocorrelation in the residuals.

FIGURE 5

FIGURE 5. (A) Scatterplot of observed and modeled (Eq. 1) DSG at Parma OBS for the calibration sub-set (1938–1990), with regression line (red line), the bounds showing 90% confidence limits (pink colored band); (B) histogram of residuals; (C) Q-Q plot (sample vs. theoretical quantile values).

The model residuals approximate a normal distribution (Figure 5B; normality test, p > 0.05, Jarque and Bera, 1981) and the Q-Q plot (Figure 5C) exhibits a distribution of sample-quantiles around the theoretical line, indicating only a few biased high DSG values.

At the validation stage, the value of the R² statistic indicates that 56% of the total variability in the observation is explained by the model, while MAE is equal to 4.4 days year⁻¹, which is lower than the standard error of the estimates (5.7 days year⁻¹), while the Kling-Gupta index is 0.6 (equal to the value obtained with the calibration set). Also in this case, the Durbin-Watson statistic (DW = 1.60, p = 0.13) indicates that there is no significant autocorrelation in the residuals. The timeline of DSG validation period shows the coevolution of observations and estimates (Figure 6A). The related distribution of model residuals (Figure 6A1) does not deviate significantly from normality (p > 0.05). In Figure 6B, the model estimates are compared to the decreasing trend (from 1950 to 2005) of the index obtained by De Bellis et al. (2010) for the Emilian Apennine, by averaging the snow-cover values of different mountain stations normalized with respect to the variability. We evaluated the relative performance of the DSG(Parma) model, without comparing the absolute estimates. Coevolution between the SCNI (snow-cover normalized index) and the modeled DSG illustrates a substantial agreement (r = 0.74).

FIGURE 6

FIGURE 6. (A) Timeline evolution of observed (blue curve) and modeled (Eq. 1, orange curve) DSG for the validation sub-set (1991–2018); (A1) histogram of residuals; (B) timeline coevolution of observed Snow-Cover Normalized Index (SCNI) for the Emilian Apennine mountain range (blue curve) and modeled DSG at Parma OBS (orange curve) for the period 1950–2004.

We could get a slightly better performance by replacing P_w in Eq. 1 with the square root of snow days per year during both the calibration (R² = 0.76, KGE = 0.81, MAE = 8.6 days year⁻¹) and validation (R² = 0.61, KGE = 0.52, MAE = 4.6 days year⁻¹) periods. We concluded that this improvement is not such as to justify the limitation of the reconstruction to the SDY observation period alone (1777–2018). Equation 1 thus supports a broader reconstruction perspective. In determining whether the DSG (Parma) model can be simplified, we have fitted a multiple linear regression model to describe the relationship between DSG and the three independent variables P_w, T_w and VC in Eq. 1. The analysis of variance of the model gave p < 0.05 for each independent variable, the highest p-value being 0.04 for the term VC. This means that the model cannot be simplified further and its results correspond to criteria of stability, interpretability and usefulness (after Royston and Sauerbrei, 2008).

Time-Series Reconstruction of Snow-Cover Persistence on the Ground (1681–2018 CE)

Equation 1 was used to reconstruct the evolution of DSG over the period 1681–2018 CE (Figure 7A). The reconstructed time-series was then analyzed to find out possible climate patterns explaining variation in long-term trends of DSG, and to compare contemporary with historical DSG anomalies.

FIGURE 7

FIGURE 7. Overview of days with snow on the ground (DSG) patterns at Parma OBS along the period 1681–2018 CE: (A) Timeline evolution of DSG data reconstructed by Eq. 1 (light blue curve), with over-imposed 11-years Gaussian-filtered series for the whole period (bold blue curve) and the observational period 1938–2018 CE (red curve); (B) Mann-Whitney-Pettitt test statistic with the change point of 1897 (vertical red arrow); (C) time series of the Atlantic Multidecadal Variability (AMV; Wang et al., 2017) together with a polynomial (third-order) interpolation curve; (D) wavelet-coherence spectrum of the standardized DSG and AMV time series; bounded colors identify the 0.05 significance level areas; the bell-shaped, black contour marks the limit between the reliable region and the region below the contour where the edge effects occur (a.k.a. cone of influence); black arrows show the relative phase relationship, with in-phase pointing right, anti-phase pointing left, and AMV leading DSG by 90° pointing straight down.

Overall, the temporal evolution of annual values presents a downward trend (Mann-Kendall M-K trend test p < 0.01; Kendall, 1975), which appears more prominent in recent times (Figure 7A, blue curve and related Gaussian fitting; Figure 7B, Mann-Whitney-Pettitt statistic). We have used different statistical methods to identify distinct climatic patterns in the long-term DSG series. In fact, different test statistics are variously sensitive to change points located at the beginning, in the middle, or at the end of a time series (Martínez et al., 2009; Toreti et al., 2011), and a combination of statistical methods is considered to be most effective to track down change points (e.g., Wijngaard et al., 2003). The application of test statistics by Pettitt (1979) and Buishand (1982) suggests the existence of a change point in the year 1897. Since the late 19th century, the Atlantic Multidecadal Variability (AMV) has experienced a significant upward trend (Si et al., 2020), coinciding with the onset of the warming period (Figure 7C). Based on terrestrial proxy records from the circum-North Atlantic region, the AMV reconstruction by Wang et al. (2017) exhibits pronounced variability on multidecadal time-scales. This multi-decadal climate mode originates dynamically in the North Atlantic Ocean and propagates throughout the Northern Hemisphere via a suite of atmospheric and oceanic processes (Wyatt et al., 2012; Wang et al., 2018b). Sutton and Dong (2012) argued for the existence of a causal link between a positive phase of the AMV and drier conditions over the Mediterranean Basin. It is also known that winter precipitation in northern Italy is mainly caused by large-scale fronts of North Atlantic and Mediterranean synoptic low-pressure systems, which produce moderate but continuous precipitation (Hawcroft et al., 2012). Interannual to multi-decadal variability of this kind of precipitation is mainly modulated by the North Atlantic Oscillation (NAO; Pinto and Raible, 2012; Gómara et al., 2014; Gómara et al., 2016), describing the fluctuations in the difference of atmospheric pressure at sea level between the Icelandic Low and the Azores High.

Other test statistics detected a change point in 1830 (SNHT-double shift; Alexandersson and Moberg, 1997) or 1967 (Worsley likelihood; Worsley, 1986; penalized maximal t-test; Wang et al., 2007; SNHT-single series, Alexandersson, 1986), which merely support the idea of a long transition period going from the final phase of the Little Ice Age (LIA; ∼1300–1850 CE; Miller et al., 2012) to the most recent warming. These different statistically-relevant years provide a loose picture of climate-related DSG variations with changing climate patterns, where the cold conditions of the LIA are still dominating after the end of the Dalton minimum of reduced solar activity (∼1790–1830; Wagner and Zorita, 2005) until toward the end of the 19th century, but in the process of evolving into an incipient warming that becomes noticeable later in the 20th century (mostly by about 1960s). We refer hereafter to the year 1897 as a relevant change-point year, which is also the starting point of more erratic weather conditions. The long-term mean value of DSG (Table 1) is equal to 28 days year⁻¹ (±16.8 days year⁻¹ standard deviation) until the change-point year (Figure 7B) and to 21 days year⁻¹ (±14.8 days year⁻¹) after that point for the continuous predicted time series (Figure 7A, blue curve), while it is equal to 20 days year⁻¹ (±11.9 days year⁻¹) for the joined modeled and observed series (Figure 7A, joined blue and ochre curves after the year 1897).

TABLE 1

TABLE 1. Descriptive statistics of the DSG (days with snow on the ground, day year⁻¹) time series for two climatic periods.

Considering that these two series, as well as the predicted and observed series for the period covered by observations (1938–2018 CE), are not statistically different (paired Student-t test, p-values >0.05), our analysis was based on the modeled values only. Table 1 also shows that the median and the 75th percentile values were higher over the period 1681–1897 CE, while more extreme values (95th percentile) were higher after the change point. The trend toward anomalously warm conditions over the 20th century and the beginning of the 21st century indicates that increasing temperatures represent the main factor triggering the decline of DSG, in particular over the most recent five decades where DSG have been subject to further decrease. This is in agreement with the contraction of snow-cover extent recently observed across Siberia and North America (e.g., Musselman et al., 2017; Zhang and Ma, 2018). It is striking that snow covered the ground over 145 days in the year 1830 and only over 99 days in the year 1709 (Figure 7A). The latter proved to be the coldest in memory, with a particularly long and hard freeze in a winter passed to the historical chronicles because of its cold and snowstorms. In the central Mediterranean, the severity of this winter was accompanied by ice and frozen water bodies like lakes and rivers (Glaser et al., 1994; Grove and Conterio, 1994), while snowfall continued from late December to mid-February, after which rain and snowmelt caused several rivers to overflow (Diodato et al., 2019). However, the 1829–1830 winter was not only described as harsh in all parts of Europe, but also early, with long and heavy snowfall from mid-November to the end of March (Corradi, 1865–1890, vol. III, p. 403). In January 1830, even the press was suspended because the roads, due to the abundance of snow, had not allowed the arrival of newspapers. The year 1830 was rather the end of a period that extended over important parts of the LIA, during which Europe experienced predominant cooling (Xoplaki et al., 2001). In fact, winter 1835–1836 was also remarkable for Milan and its region but not quite as severe as others (Stella, 1836, pp. 350–351):

The year 1684, at the beginning of our series, is also distinguished by a high number of days with snow on the ground (88 days). In this regard, Corradi Annals’ (1865–1890, vol. II, p. 257) define 1684 as the year in which snow covered the ground until after Easter:

After Xoplaki et al. (2001), the severe conditions during the winter 1683–1684 can be attributed to the presence of an extended high-pressure system over central and western Europe, with its center over north Iberia. This distribution of atmospheric pressure centers in Europe led to persistent northwesterly and/or northeasterly circulation over the Balkans and the eastern Mediterranean. These cold air outbreaks caused low or extremely low temperature conditions connected with precipitation events (rainfall and/or snowfall). Also noteworthy are the years 1947 and 1963, with snow on the ground over 79 and 73 days, respectively. January 1947 was recognized as historical for Italy, with very strong waves of frost and snow all over the country. The coldest winter in the post-World War II period was, however, 1962–1963: started on December, it continued with an escalation of continuous Siberian cold spell, alternating with Atlantic phases, until the month of March, which was equally rigid. In Britain, although the Thames had not frozen as in previous centuries, that winter turned out to be not only the coldest of the century but among the coldest ever. After that winter, the DSG median dropped further (only 15 days per year), and with it also the 95th percentile reached the lowest values ever (30 days per year). These trends are confirmed by the decrease in snow-cover extent in spring, as observed in Eurasia. Although the extent of the autumn snowpack has been limited in recent decades, and the winter trend has remained unchanged in Europe and Asia, a decrease in the snowpack extent and a greater variability in the transition seasons (spring and autumn) have been documented at the hemispheric scale since the 1920s (Krasting et al., 2013). Such snow-cover reduction was particularly significant in the mid-latitudes (40°–60°) of the Northern Hemisphere (Brown, 2000). Accordingly, Déry and Brown (2007) noted that the mean monthly snow-cover extent strongly anti-correlates with the mean air temperature in the Northern Hemisphere. This occurs particularly in April–June, when extended snow-cover duration is accompanied by intense incoming solar radiation.

We used the wavelet-coherence analysis for highlighting the time and frequency intervals when DSG and AMV have a significant interaction (Grinsted et al., 2004). The wavelet coherence method offers a bivariate extension of the wavelet analysis to identify regions with large common power in the time-frequency domain of the two time series, and reveals information about their phase relationship (which may be suggestive of dependency; Maraun and Kurths, 2004). The wavelet-coherence spectrum of the two time series (Figure 7D) displays sporadic high-frequency periodicities around 11 years. The quasi-11-year sunspot cycle is a main feature of solar activity (as identified by Wolf, 1852; Wolf, 1853). However, the 5% significance level is a not a reliable indication of dependency for erratic periods in the long-term behavior although Italian Alps snow-cover trends have been observed to oscillate with 11-year periods (Valt and Cianfarra, 2010). The region of the spectrum at lower-frequency periodicity of ∼60 years is quite extensive and suggestive of a causal AMV-DSG relationship emerging at the onset of the warm period (across the change-point year of 1897 with AMV leading in time), which appears to be linked to North Atlantic internal ocean-atmosphere variability (e.g., Knudsen et al., 2011; Mazzarella and Scafetta, 2012). However, the significant area falls outside of the cone of influence where edge effects become important and prevents this analysis from obtaining a robust interpretation in this sense.

Influence of Atmospheric Circulation Patterns of Snow-Cover Periods

The occurrence of snow events in Parma is to some extent related to the presence of a low-pressure system over the central Mediterranean area, which is a favourable configuration for several snow days and persistent snow. In order to discover if different climatic periods are somewhat the result of changes occurred in the continental atmospheric circulation, we created composite plots of sea-level pressure for the winters of the late Maunder Minimum (Figure 8A), the late LIA until the change-point year (Figure 8B) and the modern warming (Figure 8C). These climatic sub-periods reflect many dominant types of large-scale atmospheric circulation patterns, which are responsible for wintertime conditions over the Mediterranean region. During the first period, 1681–1715 (Figure 8A), the atmospheric pattern shows that the central Mediterranean area is characterized by a winter depression, which attracts cold arctic air masses from the north/northeast of Europe. In the period 1716–1897 (Figure 8B), a similar atmospheric circulation pattern continued to favor the transition of cold air masses from northern Europe though with a shallower depression center lessening snowfall and snow-cover duration. There were certainly less frequent and severe cold spells during the late LIA than in the Maunder Minimum. However, the late LIA was also characterized by some severe winters. To give an example, taken from Corradi (1865–1890, vol. II, p. 623), in 1782 the whole of Italy was besieged by a deadly bad weather:

FIGURE 8

FIGURE 8. Mean reconstructed sea-level pressure (SLP) maps over southern Europe over the winters (December to March) of three periods: (A) late Maunder Minimum (1681–1715 CE); Late Little Ice Age (1716–1897 CE); (C) Modern Warming (1897–2008). L = low-pressure system (cyclone), H = high-pressure system (anticyclone). Maps are authors’ own elaboration using Climate Explorer (https://climexp.knmi.nl/start.cgi) [Accessed September 10, 2020] on data from: (A,B)Luterbacher et al. (2002); (C) from Küttel et al. (2010), which stops at year 2008.

The exit of the LIA generally offered less inclement winter weather. The map of the third period, 1898–2008 (Figure 8C), shows the squeezing and reabsorption of the depression center by two anticyclones, which is indicative of snow weather (both mean and extreme snowfall) becoming less common. The two high-pressure systems (Azores high to the west and Russian high to the east) prevented the Mediterranean from prolonged snow events of certain importance, as compared to the previous periods.

Conclusion

The annually-resolved DSG (days with snow on the ground) time-series reconstruction documents variations at Parma OBS since 1681. With this reconstruction, unprecedented in time length, we have deepened our understanding of the characteristics of snowfall in the Po valley (northern Italy) and reported the consistency of our model-based reconstruction of DSG with historical data. That consistency suggests that the reconstructed historical signal may be representative of real multi-decadal variations. In fact, our 339-year long time series of DSG shows a particularly marked downward trend in recent decades, after the change point detected in 1897, and suggests that important shifts in mean DSG values and their variability may be resolved by extended reconstruction. Interannual and inter-decadal variations are evident in the reconstruction. While the causes of the decrease observed in DSG are manifold and varied, two factors are recognisable: an increase in seasonal temperatures (Dobrovolný et al., 2010; Serquet et al., 2013) and a contraction in the variability of snow precipitations from 1 year to the next. With respect to the latter, we highlight a general shortening of snow duration due to a delay in the start date of snow cover in autumn and an advancement of the end date in spring, associated with a retreat of glaciers, mainly due to a decrease in winter precipitation (Vincent et al., 2005; Huss et al., 2008). However, other factors besides the increase in temperature may have caused a decrease in the snow-cover duration. Already in 1881, on the occasion of the royal meeting of the Accademia dei Lincei (“Lincean Academy”) in Rome (Italy), the Italian Catholic priest and geologist Antonio Stoppani (1824–1891) in a speech Sull’attuale regresso dei ghiacciai sulle Alpi (“On the current regression of glaciers in the Alps”) provided as a main cause of the regression of the glaciers, not the variations of temperature (even if they were of common knowledge), but the reduction of snowiness (e.g., Scaramellini and Bonardi, 2001). In fact, in that period the number of snow days had been reduced to roughly a third in just about 50 years, a worrying phenomenon, especially because at that time there were no dams and artificial reservoirs, and glaciers were among the few large reserves of water (e.g., Stoppani, 1876). Again, from 1797 to 1806 the days of snow in Milan (Italy) had been 243 (i.e., 26 on average per year), but the situation had substantially changed from 1857 to 1876 with 166 snow days (i.e., eight snow days per year). Recently, Diodato et al. (2020b) provided evidence for Switzerland that this decrease is continuing in the Alpine range.

In our study, we report the results without providing conclusive evidence regarding the causes of the observed changes in the reconstructed DSG. Our model results offer robustness against anthropogenic disturbance in urban environments. The latter is relevant for the Parma territory, 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 days year⁻¹ in the 1940s to <20 days year⁻¹ in the most recent years), which appears to be sufficiently captured by the modeled time series. We have corroborated our findings suggesting that change points in snowfall series in northern Italy can be linked to large-scale changes in the modes of climate (e.g., the Atlantic Multidecadal Variability), identifying in this way the need for future research on the subject. Though the long-term trend of the DSG series may be represented by a combination of intermittent periodicities, due to the brevity of the time interval the presence of these periodicities and their relationship with snowfall mechanisms in northern Italy should be regarded as tentative, and in need of confirmation by additional studies on extended spatial scales.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

ND conceived the study, performed the analysis and drafted the manuscript with GB, who wrote the final manuscript. CB prepared and organized the snowfall data used in the study, and contributed to scientific discussion of the article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the NTNU publishing fund support to cover article-processing charges. The authors would like to thank Iñigo Gómara (Universidad Complutense de Madrid, Spain) for his support on wavelet-coherence analysis (wavelet software was provided by A. Grinsted).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.561148/full#supplementary-material

References

Adamson, G. C. D. (2015). Private diaries as information sources in climate research. WIREs Clim. Change. 6, 599–611. doi: 10.1002/wcc.365