Recent waning snowpack in the Alps is unprecedented in the last six centuries

Snow cover in high-latitude and high-altitude regions has strong effects on the Earth’s climate, environmental processes and socio-economic activities. Over the last 50 years, the Alps experienced a 5.6% reduction per decade in snow cover duration, which already affects a region where economy and culture revolve, to a large extent, around winter. Here we present evidence from 572 ring-width series extracted from a prostrate shrub (Juniperus communis L.) growing at high elevation in the Val Ventina, Italy. These ring-width records show that the duration of current snowpack cover is 36 days shorter than the long-term mean, a decline that is unprecedented over the last six centuries. These findings highlight the urgent need to develop adaptation strategies for some of the most sensitive environmental and socio-economic sectors in this region. Snow is an important component of the environment and climate of mountain regions, but providing a long-term historical context for recent changes is challenging. Here, the authors use ring-width data from shrubs to show that recent snow loss in the central Alps is unprecedented over the last 600 years.

Snow cover dominates hydrological cycles and climate of high-elevation and high-latitude regions, representing the key interface between the atmosphere and the ground. Snow impacts the surface energy balance by altering albedo and emissivity and, together with its thermal insulation properties and meltwater input, can substantially affect glaciers and permafrost with a major influence on Earth's climate 1 . Snowpack also acts as a surface water storage reservoir 2 , driving the timing of runoff that sustains downstream environmental and human water demands. Billions of people around the world depend on these resources and mountain ranges are nowadays recognized as the 'world's water towers', providing, with glaciers and snowmelt, a substantial water supply to downstream areas during the dry and hot seasons 3 . The Alps are the most important water-supplying mountain range in Europe, with the Danube, Rhine, Rhône and Po being the main drainage basins and the Po river basin, on the southern side of the Alps, deemed the second most sensitive area, after the Rhône basin, when considering the intersection between water supply and demand for agricultural, industrial and domestic purposes 3 . These mountain environments are recognized as being at risk from the current climate crisis and will probably face dramatic changes in the future with cascading impacts not just in the high-elevation areas but also to the downstream regions 4 . Even just focusing on the snowpack, there is a range of environmental processes and socio-economic activities directly connected with the winter snow cover duration; some snow-dependent or cold-adapted species are facing a decline in abundance or a reduction in reproductive fitness 4,5 , while winter tourism and related recreation activities have been negatively impacted by the space and time reduction of the snow cover 4 .
In the Alps seasonal (November-May) mean snow depth, recorded in several hundreds of in-situ snow depth observations, experienced an 8.4% decline per decade between 1971 and 2019, with a parallel reduction of 5.6% per decade in snow cover duration 6 . To thoroughly understand how recent snow cover dynamics and associated processes are unusual and how they will be influenced by ongoing and future climate change, it is crucial to have information about long-term reliable observations or proxy series of snowpack extent and duration. The Alps are a mountain region with one of the longest traditions in weather record collection. Here, there are several instrumental pressure, temperature and precipitation series extending back to the mid-eighteenth century 7 . These valuable instrumental records, coupled with the abundance of long-lived tree species growing at high elevation, led to the development of several yearly resolved tree-ring climate reconstructions calibrated and independently verified to instrumental target data over exceptionally long periods. These reconstructions, with few exceptions [8][9][10] , focussed on summer temperature [11][12][13][14][15] , mirroring the key limiting factor and the temporal sensitivity window for tree Article https://doi.org/10.1038/s41558-022-01575-3 This attribute, coupled to its high longevity, turns ring width of common juniper into a potential proxy for snow cover duration across the Alps 26 .
To derive a climate proxy from ring width of prostrate juniper, we conducted several preliminary tests to assess the robustness of the long-term trends (Extended Data Figs. 1 and 2) but also the role of several growth peculiarities (for example, missing and wedging rings or eccentric development; Fig. 1) commonly present in this species (Methods). In parallel, we developed and then validated a method to estimate the seasonal evolution of snowpack over a century-long timescale starting from daily instrumental precipitation and temperature series (Methods) (Extended Data Figs. [3][4][5][6]. Snowpack duration was estimated during the snowy season considering two key components: a positive term related to the amount of solid precipitation and a negative one due to the melting process and related ablation of the surface snowpack mostly linked to temperature 27 . The snow degree-day model we applied is a simple snow cover model 28 ; however, simulating snow mass evolution solely from daily precipitation and temperature series permitted us to obtain a long record of snow cover duration (1834-2018), impossible to reach with more sophisticated approaches. This is crucial for a sound and reliable calibration and verification of the reconstruction model using the ring-width proxy we introduced. growth in heat-limited environments 16 . However, local measurements of snowpack often extend back only a few decades 6 , while the long resting period during the cold season joined with the mostly negligible moisture-limiting conditions of the Alps, prevents adopting tree rings as an effective proxy to reconstruct snow conditions 17 . As short-term series of observations are not adequate to describe snowpack dynamics and the full amplitude of mountain wintertime conditions 18 , the absence of any continuous long-term (centennial or more) instrumental or reconstructed record of snowpack duration for the Alps so far, prevents it being correctly placed in the right context and assessing the true nature of the present-day shrinking 6,[19][20][21][22][23][24] .

Reconstructing snow cover duration
We here address these issues by developing a 600 yr snowpack duration reconstruction based on 572 ring-width individual series from the southern Alps. The record was developed over 5 yr using living and relic common juniper (Juniperus communis L.) shrubs from a high-elevation (>2,000 m) site (Methods; Supplementary Table 1).
Unlike trees, that thanks to their typical erect growth are usually tightly coupled to prevailing free atmospheric conditions 25 , prostrate shrubs lack an erect life form. This habit, usual in junipers growing at high elevation (Fig. 1a), prevents most of the physiological activities and the growing onset as long as the shrub is beneath the snow cover. Calibration/verification of ring width with a modelled snowpack duration series is remarkably strong and stable (medians of bootstrap calibration/verification statistics are −0.688 and −0.687, respectively; Supplementary Fig. 1). The final snow cover duration ring-width reconstruction (SALP), calibrated against modelled yearly (October-September) snowpack duration (r 1834-2018 = 0.687), spans the 1400-2018 ce period and features significant temporal stability (Supplementary Table 2).
SALP shows a high year-to-year variability from late-lying snowpack, lasting almost all year round (as in 1431, 1541 or 1705), to years of <200 d of snow cover duration, that is around 2 months less than the long-term mean value of 251 d (for example, 1532, 1875 or 2012) (Fig. 2). We also identified a succession of episodes with long and short snow cover duration, including the 1440-1460 and 1780-1800 periods that represent the longest-lying 20 yr phases of reconstructed snowpack or the low-snowy decades 1940-1960 (Extended Data Fig. 7). The shrinkage of snow cover duration, starting around the end of the nineteenth century (as highlighted by the piecewise regression fit, orange line in Fig. 2), emerges as a preeminent occurrence within the last 600 yr and the first two decades of this century record the lowest point with a mean snow cover duration of 215 d, 36 d less with respect to the long-term mean. The co-occurring short-and long-term variability in SALP mostly derives from a blended signal that combined different time domains: a high-frequency signal, mostly related to cool-season precipitation and a lower-frequency one, much more closely related to yearly (autumn to summer) temperature (Fig. 3).

Comparison with other sources
As of today, no study has reconstructed yearly snow cover duration in the Alps for century time scales. Thus, a comparison with other investigations or climate reconstructions is not straightforward, due to the different temporal coverage or the different targeted signals. Nevertheless, our findings align well with recent studies examining snow cover dynamics over the same region, although over much shorter periods. The shrinking trend we reconstructed in the last five decades (−5.4 ± 2.5 d decade −1 ) greatly overlapped with that derived by ref. 6 for the same sector of the Alps, with observational data and when the full snowy period is considered (−6.67 days decade −1 ).
The blended temperature and precipitation signal, encoded in juniper ring width, mirrors the complex and dynamic interplay of air temperature and precipitation that affect snow cover duration but also glacier dynamics. The 500 yr reconstructions of the cumulative length-change for seven Alpine glaciers 29 , featured a long-term change closely matching what we found in SALP where a centennial phase of stability, from the sixteenth to the mid-nineteenth centuries, follows a rapid retreat briefly interrupted by the ephemeral recovering stage between 1960 to the early 1980s (Extended Data Fig. 8).
Distilling on ref. 8 seasonal precipitation and temperature reconstructions in the European Alps, those years exceeding in turn the 10th and the 90th percentiles for autumn to spring (September-May) precipitation and spring to summer (March-August) temperature, we obtained 10 yr which are either the most dry-warm or wet-cold events. These years correspond to a precipitation anomaly >15% and a temperature anomaly >1 °C, in relation to the respective 500 yr average of the ref. 8 reconstruction. In SALP, nine out of these ten extreme years exhibit the correct sign and seven of them (1603,1947,1990,1992 as negative anomaly and 1628, 1845, 1883 as positive) show a remarkable snow cover duration anomaly of >35 d, closely corresponding in SALP to a decile departure from the mean estimated over the same 500 yr period. Comparison with documentary records shows qualitative similarities but also contradictions: out of the five severe drought events (1540, 1590, 1616, 1718 and 1719) described for central Europe 30 , 1540, 1718 and 1719 are clearly expressed as years with short snow cover duration also in SALP. Further years with relevant anomalies in SALP, such as the very short snow cover duration seasons of 1532 and 1552 (respectively, 60 and 41 d shorter than the long-term mean) or the very long-lying snow cover of 1692 (74 d longer), find consistent documentary evidences in several close Alpine areas reporting very warm winter months and early snow melting, with early vegetation onset or, in contrast, an unusually persistent snow cover 31,32 .
Winter 1916/1917, when World War I was raging in Europe and the Austro-Hungarian and Italian armies were facing each other on the crests of the southeastern Alps, is remembered for the harsh weather conditions and the exceptional amount of snow precipitation that killed thousands of soldiers-a number comparable to that caused by enemy fire 33 . In accordance, our reconstruction shows 1917 as the most persisting snow cover year of the twentieth century with 67 d beyond the long-term mean. Conversely, the signatures of several well-known harsh winters across Europe, as in 1740 or 1985, have been poorly expressed in SALP and confirm, also with available observations for 1985, the complex modulation of snowpack from temperature and precipitation where even an extreme cold spell or some exceptional snowfalls might not necessarily induce a persisting snow cover.

Discussion
The snow amount in terms of water equivalent is by far the most widely adopted parameter in snow hydrology 34,35 , with several reconstructions centred in regions where tree-ring width can act as a moisture-or energy-limited snow proxy 17,36 . The Alps, especially at high elevation, are not a moisture-limited habitat and common juniper is well-known for its frost and drought resistance 37 . This accounts for why the proxy we introduced behaves differently. Its typical prostrate growth, with a height of the order of a few decimetres, prevents juniper from discrimination between being under a 20 or a 200 cm thick snow cover. Rather, it is mostly sensitive to the time of the crown emergence from the snowpack 38 , hence the snowpack duration signal encoded in its ring-width sequences.
Snowpack duration is a matter of a complex combination of abiotic and biotic factors where the timing, magnitude and distribution of precipitation, local terrain, solar, thermal radiation and overwinter ablation, wind relocation and loading, snow features, interception of vegetation, timing and magnitude of spring warming and rain on snow  events, may differently affect the within-year snow accumulation and melting dynamics, with remarkable differences in the year-to-year snow persistence 39 . The reliance on an assorted assemblage of factors justifies the challenges in modelling snowpack dynamics, both in space and time, especially considering that the influence and relative role of some of the abovementioned factors might fluctuate from year to year. We are aware that SALP reconstruction represents a simplification of the processes controlling snowpack persistence. Yet, considering the prominent role of both precipitation and temperature in snow cover dynamics and allowing for the nature of our proxy, a prostrate shrub that is usually beneath the snow cover for the whole cold season, we can extend the validity field of SALP far beyond the typical limited range associated with the high variability of precipitation distribution and snow accumulation in a mountain terrain. In fact, juniper ring-width chronology shows remarkable reconstruction skills with both modelled (Fig. 4) and observed (Extended Data Fig. 4c) snow cover duration series across most of the southern side of the Alpine chain at similar elevation. The strong association between SALP and temperature in the low-frequency domain suggests that the recent shrinkage of snowpack duration is mostly related to the warming occurring in the last century, whereas the year-to-year variability seems to be mostly associated with the typical variability in cold-season precipitation. The absence of any significant decreasing trend in the cold-season precipitation record, together with the recent occurrence of several record-breaking events (2001, 1977 and 1983 are the first, second and fourth highest year on the 1834-2018 record for October-May precipitation; Fig. 3) reduces the perception, of local residents, tourists and authorities, of how remarkable the long-term reduction in snowpack persistence of the last decades is [40][41][42] . However, those extreme snow accumulation events did not turn into any exceptional long-lying snowpack as current warmer temperatures are always inducing a faster melting 43 .
The Alps with their snow cover represent the 'water towers for Europe' and play a key role in both natural environmental systems and a wealth of snow/water socio-economic sectors that revolve around winter snow availability. A persistent reduction in snow amount and duration will probably produce profound effects upon these systems with severe cascading impacts on human well-being. In 2022 the southern Alps experienced an exceptional snowless winter, followed by an extreme European summer drought alongside the sudden July collapse of the Marmolada glacier in the Dolomites, with its tragic death toll. Realizing that we are experiencing, together with those short-term events, a monotonic and unprecedented reduction of snowpack duration over the last six centuries, should raise public awareness of the urgent need to develop adaptation strategies and start thinking about a reform of some of the most sensitive socio-economic sectors.
Common juniper is the most widespread long-lived plant across the Northern Hemisphere 44 ; therefore, together with the array of taxa with comparable attributes (longevity, crossdatable ring parameters and prostrate growth form), what we did for the Alpine region could be realistically applied in many snow-prone regions. Information gained with SALP can help to improve state-of-the-art climate models to accurately represent snow on the ground in land surface models and properly consider snow dynamics variability at interannual to centennial time scales. This knowledge is essential to assessing whether models can correctly represent this challenging component of the cryosphere that combines modern forcing drivers with persisting natural variability 27,45 and is crucial to surface energy balance in the climate system and the hydrological cycle.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-022-01575-3.

Meteorological records and snow persistence estimation
Snow cover data are available in a systematic way for only a few decades, usually from the 1990s thanks to automatic stations and more sporadically over a longer timespan from manual stations. However, these series are too short to guarantee a sound calibration with ring-width proxy. For this reason, we developed a method to estimate historical seasonal evolution of snowpack from daily precipitation and temperature series, which are among the few variables available on a century timescale. Temperature and precipitation data come from many years of data rescue and homogenization activity, starting from the core datasets collected and presented in ref. 48 for Italy and in ref. 7 for the HISTALP project. In recent years, with follow-up data collection, most of these climate records have been recovered and homogenized at a daily resolution and the dataset extended with hundreds of station records coming from various data rescue activities by different Italian regional environmental agencies (Extended Data Fig. 3). These series were finally subjected to the same homogenization procedure (see refs. 48 and 49 for details about homogenization procedure and ref. 50 for its performances).
To derive snow cover duration, we need to estimate both the amount of solid precipitation (SP) and the snowmelt (M), then accumulate both terms at a daily step during the whole snowy season. To this end, we first used daily precipitation data, together with daily minimum and maximum temperatures, to estimate the solid fraction (in terms of water equivalent) of precipitation, consistently with ref. 26. Then, starting with mean daily temperature data, we developed and calibrated with available snow depth data, a degree-day model to estimate the snowmelt.
To estimate SP fraction, we aimed to improve the widely adopted mean daily temperature thresholds, ranging between 1 °C (refs. 51 and 52) and 2 °C (refs. 53 and 54), used to discriminate between solid and liquid precipitation. This approach implies that temperature probably exceeds that threshold during part of the day. On the contrary, by using daily maximum and minimum temperature data, we obtained a more realistic evaluation of the precipitation solid fraction, having a better representation of the temperature range over the day.
Specifically, for the ith day with a total precipitation amount P i and daily maximum and minimum temperature values Tx i and Tn i , the SP was defined as: with threshold (th) = 2 °C as in ref. 26. In the second step, snowmelt was evaluated by applying a degree-day method. This method proved to be reliable for computing the total amount of snowmelt for time lags ranging from a week to the entire snowmelt season 55,56 . It considers that a large portion of the energy required to melt snow is provided by incoming longwave radiation and sensible heat flux, both linked to air temperature. In addition, to predict snow melting for century-long time scales, the degree-day method is often favoured over the surface energy balance ones due to its sole requirement of temperature as the input variable 57 . In these cases, temperature is used to calculate the degree-day sums for a given time period, then these values, multiplied by the snow degree-day factor (SDDF), are used for snowmelt estimation 57 . Therefore, for the ith day, the amount of M is defined as: where Tm i is the mean daily temperature of the ith day, DDth the daily temperature threshold and SDDF is the snow degree-day factor.
Reference 58 , studying snowpack evolution at the Forni Glacier (Italy), suggested using DDth = 0 °C when, as in our case, the aim is to assess the amount of snow ablation. The same authors also defined SDDF as the amount of snow (in mm of water equivalent) melted per degree-day (mm °C −1 d −1 ).
where SP sum is the total snow amount (in terms of water equivalent) accumulated over the snowy season and N is the number of days necessary for melting the whole snow cover 58 . Finally, the temporal evolution of snow depth (HS) in mm of water equivalent (SWE) is defined by accumulating day by day the difference SP i − M i over the snowy season ) and the snow season length of a specific year is defined as the number of days of snow-covered ground from 1 October of the previous year to 30 September of the current year. This approach to estimate snow-season length series for any specific location was applied to the Ventina sampling site to calibrate the juniper ring-width chronology and then extended to a 30 arcsec resolution digital elevation model covering the southern side of the Alps, to verify the spatial representativeness of the chronology itself. To this end, we interpolated daily temperature (mean, minimum and maximum) and precipitation onto any specific location following ref. 59 as already applied in ref. 26.
Besides daily temperature and precipitation series, the missing parameter we still need to apply equations (1), (2) and (4), is the SDDF for each specific location. We used HS observations available across the Alps (obtained from refs. 6 and 60 and further national and regional meteorological services, such as MeteoSwiss and ARPA Lombardia, Extended Data Fig. 4a) and interpolated the empirical SDDF values onto the site of interest.
To estimate SDDF we calculated, for each station and for each snowy season, the ratio between the estimated SP and the degree-day amount as in equation (3), both cumulated over the specific snowy season period with permanently snow-covered ground, as observed in the HS station data (from the first to the last day for which HS is continuously >0).
Out of about 500 HS stations located mostly in the southern Alps, we selected a subset of 312 stations on the basis of the quality and completeness of their HS records, preferring the recent automatic period over the manual one for its completeness and reliability (Extended Data Fig. 4a).
Because the necessary meteorological variables (temperature and precipitation) were not always available at the HS stations, a site-specific synthetic daily precipitation and daily minimum and maximum temperature series were estimated with the abovementioned procedure. These synthetic series were used to estimate SP (as defined by equation (1)) instead of using the measured HS (in particular its first difference, to calculate snow gain and loss) because (1) HS observations are expressed in cm of snow and fresh snow density is usually not available to convert them into mm of water equivalent; (2) HS day-to-day variability can be significantly biased by wind-blown accumulation inducing, together with the typical rain gauges undercatch, an overestimation of SDDF and, therefore, a shortening of the snow cover duration evaluation. Therefore, the only information we used from the HS station records is the dates of the beginning and end of continuous snow cover for each year (October-September).
Next, for each HS station site, we computed the SDDF as the median (more stable with respect to the persistence of some residual spikes in the observations) of each single snow season SDDF value. Article https://doi.org/10.1038/s41558-022-01575-3 SDDF shows, as expected, a clear dependence on elevation, with higher values at higher elevation; at fixed degree-days, the higher the elevation, the more efficient the melting process due to the usually lower residual optical depth of the atmosphere and shading effect by surrounding orography. We considered this relationship with elevation to interpolate the observed SDDF values at site level. The interpolation approach is based on the estimation of a weighted linear fit of SDDF with elevation, where weights consider the similarities in the geographical settings of each station with respect to the target site in terms of horizontal and vertical distance, slope orientation, steepness and distance from the sea 61 .

Validation of the snow persistence estimation
Both SDDF and the final snow season length estimations were validated over the HS available station data falling within the analysis area domain (194 stations; Extended Data Fig. 4a,b) with a leave-one-out cross-validation: we excluded one HS station series per time and reconstructed a SDDF value representative of each station site independently from observations.
Having interpolated the SDDF value, daily precipitation and minimum and maximum temperature series on each HS observation site, we were able to estimate the daily snowpack evolution (in terms of water equivalent) over the year, using equations (1), (2) and (4). Then, with the leave-one-out approach, we compared the site-level reconstructed SWE and observed daily HS series in terms of snow-covered and snow-free conditions (contingency table shown in Extended Data Fig. 5), being that the observed HS data (in cm of snow) and the reconstructed SWE (in mm of water equivalent) were not directly comparable. However, we did not focus on single day condition but rather on the mean observed versus reconstructed length of the snow cover for the whole season. Thus, we performed a further validation computing the correlation between observed and independently reconstructed (in leave-one-out) snow cover duration series for each HS site with at least 15 yr of available data (Extended Data Fig. 4b; 135 series out of the 194 inside the working domain): mean (median) correlation value is r = 0.75 (0.79), with lower and upper quartiles ranging from 0.69 to 0.85. As an example, Extended Data Fig. 6 shows the observed versus independently reconstructed snow season length series for some of the longest HS stations series distributed across the southern Alps.

Ring-width data collection and processing
Samples were collected in the central Italian Alps, at the Ventina valley (46° 18ʹ N, 9° 46ʹ E), within or above the current treeline, between 2,100 and 2,400 m above sea level, by randomly selecting the shrubs and saw-cutting a basal disc from one of the main prostrate stems. The area is homogenous in terms of slope, aspect and soil features, with scattered trees and almost no human-related disturbances such as livestock grazing, fire, logging and so on. Soils were typically shallow and with a high gravel matrix.
At the laboratory, discs were sanded with progressively finer grit sandpaper for a clear visualization of the annual rings. In some cases, to enhance readability, we prepared the most critical samples for microscopic analysis through cutting thin sections (12-15 μm) with a microtome, staining them with safranin and permanently fixing the slides with a mounting medium 62,63 . Considering the typical irregular and lobate growth form of juniper stems at this elevation (Fig. 1b), for a more reliable representation of ring-width series, we measured two to four radii per disc. Crossdating was accomplished applying a three-step procedure: first visually comparing the different radii within each sample, second considering the mean individual series from different samples 64 and, last, checking for dating and measurement errors with the COFECHA programme 65 . Annually resolved site chronologies were obtained from the crossdated ring-width series using the ARSTAN programme that was specifically developed to remove any biologically induced age-trends and transient disturbance pulses present in raw ring-width series and to enhance the common high-to low-frequency variability probably associated to climate.
J. communis, as most of the taxa within this genus, usually features slow and irregular growth with erratic cambial activity and missing rings, not always associated with disturbances or extreme climatic events ( Fig. 1c and Supplementary Fig. 2). Yet, this species highlights a negligible age-related trend in ring width compared to erect tree species as represented by aligning the raw ring-width measurement series by cambial age (Extended Data Fig. 1). We tested different approaches to remove the age-related and disturbance trends: (1) a deterministic approach, using a generalized exponential 'Hugershoff' curve; (2) a stochastic one adopting an age-dependent spline with the initial value of 50 yr and (3) a biological approach using the regional curve standardization (RCS). With the first two methods, a specific individual detrending curve is computed for each ring-width series, whereas with the RCS a single mean growth curve is obtained, realigning by biological age, averaging and then smoothing with an age-dependent spline 66 , resulting in a growth ageing curve from all ring-width time series. We also tested an array of six different stratified samples, filtering some attributes able to potentially introduce bias in the RCS ring-width detrending and in the final mean chronology computation. These groups were created with samples with/without missing rings, with/ without the pith, with just relic or living shrubs. The growth curves (Extended Data Fig. 1) and resulting chronologies computed with the different approaches (Extended Data Fig. 2) were very similar, highlighting that these different attributes were not able to introduce any significant bias into the mean growth series. Therefore, we built the final chronology collating all the crossdated series and, to preserve at most the low-frequency signal retained in the ring-width series, we opted to apply the RCS detrending technique 67,68 . Individual series were therefore standardized by dividing observed by expected values and the final mean chronology was computed using the biweight robust mean that automatically discounts the influence of outliers in the computation of the mean, reducing the variance and bias caused by outliers 69 (Fig. 1d).
On the basis of previous investigations highlighting the key role of cold-season precipitation on common juniper ring-width variability, various calibration trials of the chronology were conducted to tune the best time window of snowpack persistence ( Supplementary Fig. 3). We then adjusted the variance and mean of the ring-width chronology against instrumental targets to avoid the typical loss of amplitude due to regression error 70 . Reconstruction skill has been assessed by considering the explained variance (R 2 ), reduction of error statistic and coefficient of efficiency statistics, together with the more stringent bootstrapped transfer function stability test 71,72 . The last evaluates, between calibration and verification periods of equal length, the stability for intercept, slope and explained variance of different subperiods over 100,000 iterations (Supplementary Table 2). The confidence intervals of the reconstruction were computed considering the standard error (±2 × root mean square error) derived from the calibration against modelled snow cover duration over the whole 1834-2018 period. Lastly, we applied the extreme value capture non-parametric test to assess the loss of skill due to scaling; it provides the statistical significance of a given number of the extreme years correctly 'captured' beyond the upper and lower 10% thresholds of the measured climate data 70 (Supplementary Fig. 4).
To assess the contribution of temperature and precipitation in the different time domains, we filtered and contrasted the frequency components from the juniper ring-width chronology and from October-September temperature (representing both the biological and hydrological year) and October-May precipitation (mostly falling as snow) series using the fast Fourier transform 47 . We applied both a Article https://doi.org/10.1038/s41558-022-01575-3 low-and a high-pass filter to block all frequency components below/ above the cut-off frequency of 10 yr (Fig. 3).

Data availability
The postprocessed ring width and modelled instrumental snow duration reconstructions that support the findings of this study and the final snow cover duration reconstruction (SALP) are available from the open repository: https://doi.org/10.5281/zenodo.7330950 ref. 73. Most of snow depth data come from the open repository 'Snow cover in the European Alps: Station observations of snow depth and depth of snowfall (v.1.1)', https://doi.org/10.5281/zenodo.4064128 ref. 60. Meteorological records (temperature, precipitation and additional snow depth) that support the findings of this study are available from the Global Surface Summary of the Day, https://www.ncei.noaa.gov/ data/global-summary-of-the-day/, the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA), http://www.scia.isprambiente.it/ wwwrootscia/Home_new.html, the ARPA Piemonte, https://www.arpa. piemonte.it/rischinaturali/accesso-ai-dati/annali_meteoidrologici/ annali-meteo-idro/banca-dati-meteorologica.html, Meteotrentino, https://www.meteotrentino.it/#!/content?menuItemDesktop=111, MeteoSwiss, upon registration on IDAWEB, https://gate.meteoswiss. ch/idaweb/login.do, the Italian Air Force, the Agricultural Research Council-Research unit for Climatology and Meteorology applied to Agriculture (CREA-CMA), the Società Meteorologica Italiana (SMI), the Agenzia Regionale per la Protezione dell'Ambiente for the regions Valle d'Aosta, Piedmont, Lombardy, Liguria, Veneto and Emilia Romagna, the Lombardy Centro di Monitoraggio Geologico (CMG), Meteotrentino, the Autonomous Province of Bolzano, MeteoSwiss, the Zentralanstalt für Meteorologie und Geodynamik (ZAMG) and the Slovenian Environment Agency but restrictions may apply to the availability of these data which were used under license for the current study and so are not always publicly available. Data are, however, available from the authors upon reasonable request and with permission of the aforementioned institutions.
Article https://doi.org/10.1038/s41558-022-01575-3 Extended Data Fig. 4 | The snow depth station network and the performance of the modelled and ring-width records. a, geographical distribution of the snow depth (HS) record station networks considered in this study; blue (white) dots are the stations within (outside) the spatial domain, delimited by the red lines in all the maps, where we carried out the analyses. b, validation of the modelled snow cover duration computed through the correlation between observed and independently modelled (with a leave-one-out validation technique) snow cover duration series for each HS station with at least 15 years of available data. Empty squares highlight the HS stations with less than 15 years of data or outside the spatial domain. c, correlation between observed snow cover duration series and juniper ring-width indices from the Ventina site. Bigger and black-bolded squares are the HS stations within 1500 and 2500 m a.s.l., all the others are outside this elevation range. As in (b) we only consider the HS sites with at least 15 years of data. The white circle in all the plots highlights the position of the study site.