Recent Increases in Winter Snowfall Provide Resilience to Very Small Glaciers in the Julian Alps, Europe

: Very small glaciers (<0.5 km 2 ) account for more than 80% of the total number of glaciers and more than 15% of the total glacier area in the European Alps. This study seeks to better understand the impact of extreme snowfall events on the resilience of very small glaciers and ice patches in the southeastern European Alps, an area with the highest mean annual precipitation in the entire Alpine chain. Mean annual precipitation here is up to 3300 mm water equivalent, and the winter snow accumulation is approximately 6.80 m at 1800 m asl averaged over the period 1979–2018. As a consequence, very small glaciers and ice/ﬁrn patches are still present in this area at rather low altitudes (1830–2340 m). We performed repeated geodetic mass balance measurements on 14 ice bodies during the period 2006–2018 and the results show an increase greater than 10% increase in ice volume over this period. This is in accordance with several extreme winter snow accumulations in the 2000s, promoting a positive mass balance in the following years. The long-term evolution of these very small glaciers and ice bodies matches well with changes in mean temperature of the ablation season linked to variability of Atlantic Multidecadal Oscillation. Nevertheless, the recent behaviour of such residual ice masses in this area where orographic precipitation represents an important component of weather ampliﬁcation is somehow different to most of the Alps. We analysed synoptic meteorological conditions leading to the exceptional snowy winters in the 2000s, which appear to be related to the inﬂuence and modiﬁcation of atmospheric planetary waves and Arctic Ampliﬁcation, with further positive feedbacks due to change in local sea surface temperature and its interactions with low level ﬂows and the orography. Although further summer warming is expected in the next decades, we conclude that modiﬁcation of storm tracks and more frequent occurrence of extreme snowfall events during winter are crucial in ensuring the resilience of small glacial remnants in maritime alpine sectors.


Introduction
Ice bodies in the Julian Alps are typically in the smallest size category of glaciers, with several glacial and nival ice patches [1] and one small mountain glacier [2]. These very small ice masses have gained an increasing scientific importance in the last few decades e.g., in the studies of Mediterranean [3] glaciation during the Pleistocene and the Holocene [4] and elsewhere. Small glaciers represent relevant component of the cryosphere contributing to landscape formation, local hydrology and sea level-rise [3]. At a regional scale like the European Alps, residual ice masses may occupy a significant volume fraction and their omission could result in errors of ±10%. For this reason, several authors stressed to the production of accurate inventories considering all ice bodies down to 0.01 km 2 or smaller in order to reduce uncertainties [5,6]. Despite their small size, they are widespread and likely account for about 80-90% of the total number of mid-to low-latitude ice bodies in mountain ranges. For this reason, it is important to understand their behaviour [7].
The Alpine landscape is rapidly changing due to global warming and larger ice masses continue to disintegrate into smaller ice bodies (e.g., [2,[8][9][10]). In the future, alpine glaciers will be likely characterized by a larger number of these very small ice bodies, especially in karstic terrains. Indeed, the cryosphere is rapidly declining in most, if not all, high elevation environments on Earth [11] and anthropogenic climate change is thought to be the main cause in accelerating this process over the next few decades. In the European Alps a rapid glacier retreat and downwasting has been observed in the early 21st century with regionally variable ice thickness changes of −0.5 to −0.9ma −1 and an estimated 2000-2014 mass loss of 1.3 ± 0.2 Gt a −1 [12]. The Equilibrium Line Altitude (ELA) in the Alps rose by 170-250 m in the last few decades [13], and in the Julian Alps, the Environmental ELA (envELA; [14]), the temperature-precipitation-dependent ELA [15,16], rose from 2350 to 2750 m between the early 1980s to the beginning of the 2000s [17]. The main cause for glacier retreat in the Alps is related to longer and warmer summers modulated by winter precipitation anomalies [11,13,[18][19][20]). Nevertheless, the magnitude of these impacts is both site-and catchment-specific and depends on the resilience and buffering capacities of the system, including the area occupied by glaciers and permafrost, basin hypsometry, groundwater systems, vegetation type and cover, and ecosystem types [21][22][23]. Such ongoing major and unprecedented glacier decline strongly affects also the integrity and value of several World Heritage sites [20].
Very small glaciers are very sensitive to climatic change and subsequently their areas may change significantly on a decadal, and even on a yearly basis [24]. Nevertheless, small glaciers and ice patches in some sectors of the Alps as well as in some Mediterranean areas have surprisingly shown a slowdown in the reduction rate from the mid-2000s even though summer temperatures continued to rise at double the rate of the global average [2,17,[25][26][27][28][29][30]. This is especially true in more maritime environments where mean annual precipitation (MAP) is very high or where the topography favours snow accumulation [31][32][33] Based on high-resolution and high quality daily precipitation measurements, the Julian Alps is considered as a mesoscale hot spot of heavy precipitation [34] where the highest alpine MAP is recorded [35]. Furthermore, total precipitation has different origins over the northern and southern flanks of the chain: precipitation is more frequent but with lower intensity to the north, while much of the total precipitation to the south is due to rare intense weather events [34]. Local feedbacks and synoptic conditions can favour or even moderate such events.
Variation in the large-scale atmospheric circulation, in part associated with climatic modes of variability such as the North Atlantic Oscillation (NAO) and the Atlantic Multidecadal Oscillation (AMO), produces different synoptic conditions (Northern Hemisphere teleconnection patterns, cyclones) leading to a high interannual variability of temperature and winter precipitation over the Alps [36,37].
NAO is an irregular fluctuation of surface atmospheric pressure over the North Atlantic Ocean involving the subtropical high located over the Azores islands and the semipermanent subpolar low pressure system centered over Iceland [38]. In wintertime the NAO can be considered as the most energetic mode of climate variability in the northern hemisphere that can also manifest itself in extreme weather events [39].
AMO explains multidecadal sea surface temperature (SST) variability in the North Atlantic Ocean from 0 • to 70 • N, with a cycle of 65-70 years and 0.4 • C range between extremes [40]. Accordingly, as one of the most prominent modes of climate variability, AMO also influences Europe's climate [41]. For instance, it has already proven that AMO has impact on severe cold events across Europe as the positive phase of the AMO manifests in more frequent negative NAO, thus more blocking events associated with extreme weather conditions might occur [42]. The primary influence of AMO on alpine climate is on mean air temperature, while its influence on precipitation is less clear [43,44].
Precipitation, both in winter and summer, is strongly influenced by the position of blocking high [45]. Long-lasting high pressure system, which blocks the moisture air or Atlantic winter cyclones penetrating from the west, leads to stable conditions along with the downwelling motions [46]. Accordingly, during wintertime within a region where weather conditions are dominated by an anticyclone, over the basins surrounded by mountains, sharp temperature inversion can occur [47,48]. Furthermore, It has also been recognized already in the late 19th century that specific cyclone tracks can lead to abundant large-scale precipitation and even snowstorms in the region of central Europe [49]. It is also interesting to note that van Bebber's type cyclones (Vb, fifth subcategory of cyclone that follows van Bebber's track, [49]), which move from the east-southeast of Europe into the inner part of Europe can be related to heavy precipitation (95th percentile) events. These relatively rare events can be associated with heavy precipitation manifested over the northern slopes of the Alps.
The aims of this paper are (1) to report mass balance observations over the period 2006-2018 for very small glaciers in the Italian part of the Julian Alps (southeastern European Alps) obtained from multiannual geodetic surveys performed by using LiDAR techniques, theodolite measurements and GIS analysis; (2) to link the change in mass balance with associated weather patterns and their possible local amplifications, with a focus on extreme snowy winters observed over the same time span.

Study Area
The Julian Alps, designated in 2019 from UNESCO (United Nations Educational, Scientific and Cultural Organization) as a Biosphere Reserve (https://en.unesco.org/biosphere/ eu-na/julian-alps-italy; last accessed 15 July 2020), occupy 1853 km 2 , spanning west to east across the Italian-Slovenian border ( Figure 1). The mountain range is characterized by Upper Triassic carbonate succession of dolostone (Dolomia Principale) and limestone (Dachstein limestone), and glaciokarst landscape [50] is widespread. The drainage divide between the Adriatic Sea and the Black Sea crosses Julian Alps in the range of Mount Canin-Kanin (2587 m asl; Figure 1), which represents an important groundwater storage in the area [51,52]. Mount Montasio (2754 m asl) is the highest peak in the Italian side of the Julian Alps. The area is recognized as the wettest in the Alps with quite constant precipitation contributions from all seasons [35]. The 1981-2010 climatology of the Mt. Canin-Kanin area has been reconstructed by [28] at an elevation of 2200 m asl (AWS Canin, Figure 1a). MAP is estimated at 3335 mm water equivalent (w.e.) with February the driest month (180 mm w.e.) and November the wettest (460 mm w.e.). Mean annual air temperature (MAAT) is 1.1 ± 0.6 • C. February is the coldest month (−6.0 ± 2.9 • C) and July the warmest (9.2 ± 1.3 • C).
Average winter snow accumulation of 6.80 m was measured between December and April from 1972-2018 at 1830 m a.s.l. at the Gilberti nivological station (Figure 1a). Snow thickness on the ground is also recorded at Livinal Lunc weathter station (Figure 1a). The present mean envELA is estimated at about 2700 m asl [17], which means it is at present located above the highest peaks except for Mount Montasio. From the end of the little ice age (LIA; about AD 1865), glaciers in the Julian Alps have retreated significantly, and the majority changed in status from glaciers to either small glaciers, very small glaciers, ice patches of glacial origin [1], buried ice patches (permafrost patches), or simply no ice at all [17,31]. Since the LIA peak in ice extent, the volume decreased by about 96%. In 2012, 19 very small ice bodies still existed in the study area, having areas from 2·10 −3 km 2 to 71·10 −3 km 2 and estimated volumes from 0.002·10 −3 km 3 to 0.822·10 −3 km 3 [31]. The only ice body retaining the characteristics of a mountain glacier is the Montasio West, while the Average winter snow accumulation of 6.80 m was measured between December and April from 1972-2018 at 1830 m a.s.l. at the Gilberti nivological station (Figure 1a). Snow thickness on the ground is also recorded at Livinal Lunc weathter station (Figure 1a). The present mean envELA is estimated at about 2700 m asl [17], which means it is at present located above the highest peaks except for Mount Montasio. From the end of the little ice age (LIA; about AD 1865), glaciers in the Julian Alps have retreated significantly, and the majority changed in status from glaciers to either small glaciers, very small glaciers, ice patches of glacial origin [1], buried ice patches (permafrost patches), or simply no ice at

Topography-Related Glacier Data and Analysis
We used different geospatial data, such as LiDAR, theodolite-generated ground control points (GCPs) and digital orthophotos dated to the period 2006-2018 in order to derive outlines of all ice bodies in the Italian sector of the Julian Alps (n = 33) and calculate their surface elevation change and volume change during this period. We also calculated annual and 12-year cumulative geodetic mass balance of the Canin East GIP for the period 2006-2018. This GIP was chosen as a case study for the annual mass balance measurements because there is a long historic record of its front variations [17] and recently acquired internal structure and densities of different units within this ice body [53]. We used ESRI ArcGIS ® and the WGS1984 UTM Zone 33N spatial reference system for mapping and analyzing the geospatial data.

LiDAR Data Acquisition and Processing
Airborne LiDAR (Light Detection and Ranging) data were acquired over the study area in the years 2006, 2011, 2013, 2015, 2016 and 2018, with an average point density greater than 4 points per square meter. All the surveys were carried out in the same period of the year, between 13 September and 1 October. Technical specifications of the data acquisition process are reported in Table 1. The point clouds were then exploited to create Digital Elevation Models (DEMs) with a resolution of 1 m. In particular, the software employed for this task, TerraScan by TerraSolid, uses LiDAR points to derive a Triangulated Irregular Network (TIN) through Delaunay triangulation that is subsequently resampled onto a grid of specified resolution.

Theodolite-Generated Ground Control Points (GCPs)
Manual field surveys were performed using a total station instead of GPS due to no or patchy network real time kinematic (NRTK) available correction and the difficulty of obtaining a sufficient number of satellites. To do this, five ground control points (GCPs) were previously located in the vicinity of the GIP to assure repeatable measurements with a total station. The total station was a Stonex R9 DR1000 robotic. One of the GCPs is located just on the Canin automatic weather station (AWS) that ensures good visibility even in case of large amount of snow on the ground. Points on the GIP were surveyed using a pole-mounted 360 • prism. The furthest surveyed point is about 300 m from the GCP, and with this configuration, the accuracy of the survey is centimetric.
GCPs (Table 2) were exploited to create digital elevation models (DEMs) with a resolution of 1 m for every year between 2012 and 2018. Different interpolation methods, such as natural neighbour, spline, inverse distance weighted and topo to raster, were tested to derive DEMs from GCPs. We computed a DEM of Difference (DoD) by subtracting the LiDAR derived DEM from GCPs-derived DEM for the years 2013, 2015, 2016 and 2018, when both LiDAR and theodolite generated GCPs data were available. Average elevation difference and standard deviation were computed for every DoD. Both natural neighbour and spline interpolation methods gave on average the smallest average elevation difference (0.26 m) between the LiDAR derived and GCPs-derived DEMs (Table 3). Moreover, natural neighbor interpolation method gave on average the lowest standard deviation (0.36 m). Therefore, GCPs-derived DEMs for years 2012, 2014 and 2017, for which there is no LiDAR data available, were generated using natural neighbour interpolation method and used for further GIP-related computations. For years when both theodolite-generated GCPs and LiDAR data are available, we favour the LiDAR derived DEMs for GIP surface elevation change calculations due to higher point density and consequently higher accuracy of DEM. GCPs served also for quantifying ice-body dynamics for the period 2012-2017.

Ice Body Measurements
We manually digitised the surface area of all the glaciers and ice patches in the Italian Julian Alps for the years 2006 and 2018 with the use of digital orthophotos and 1 m resolution LiDAR-derived shaded relief. The associated absolute error by using this method, according to Fitzharris et al. (1997) [54] is about ±25 m 2 . Abermann et al. 2010 for glaciers smaller than 1 km 2 (which is 2 orders of magnitude larger than our ice bodies) indicate a larger error of ±5%. The surface area of the Canin East GIP was mapped for the year 2006 and for every year between 2011 and 2018 using digital orthophotos, 1 m resolution LiDAR derived shaded relief or theodolite-generated GCPs taken at the GIP margin ( Figure 2). Each DEM was clipped by the corresponding GIP area so that all cells in the DEM represent snow-, firn-, and ice-covered areas. We then computed a DoD of sequential DEMs by subtracting DEM at a later time from a DEM at an earlier time in order to obtain the surface elevation change of the Canin East GIP for every year between 2011-2018, a 5-year Average surface elevation change of each DoD was then multiplied with the corresponding surface area, which yielded a volume change.
Atmosphere 2021, 12, x FOR PEER REVIEW 8 of 26 using a floating-date systems [56], where the surveys were carried out close to the end of the ablation period (see Tables 1 and 2).

ERA Interim Dataset
The meteorological variables (wind fields at 10 m above surface, sea level pressure, precipitation) used in this study in order to analyse synoptic conditions are provided by the ERA-Interim dataset [57]. The ERA-Interim is an atmospheric reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) for the whole globe. The data assimilation system used to produce ERA-Interim includes a 4dimensional variational analysis (4D-Var) technique with a 12-h analysis window. Additionally, the ERA-Interim dataset is freely available for the period 01.01.1979-31.08.2019 (noting that the ERA-Interim dataset has been superseded by the ERA5 dataset, [58]) at approximately 80 km grid resolution and on 60 vertical levels from surface to 0.1 hPa, which makes this dataset ideal for analysing synoptic situations leading to abundant precipitation or long-lasting rainless events over a specific region in Europe. Figure S1 reports Unlike the volume change calculations, where the yearly area change was fully taken into account, we used a common area to compute the geodetic mass balance of the Canin East GIP. The common area is the intersect area from all observational years in the period 2006-2018 and, moreover, does not include the uppermost part of the GIP that is mainly debris-covered ( Figure 2). Thus, the common area is the most stable area of the GIP and therefore the most appropriate for the mass balance calculations. Geodetic mass balance was computed by multiplying average surface elevation change with the density of the GIP material. The density of this ice body was previously estimated to 791 kg m −3 by [53] using ground penetrating radar in 2011, and this value was applied to compute the geodetic mass balance. Due to larger quantities of residual winter snow/firn in years 2013-2015 and 2018, a density of 650 kg m −3 was used in mass balance calculations for these years [53][54][55] were also obtained. The latter was calculated also for all other ice bodies in the Italian Julian Alps using the intersect area, and the conversion factor of 791 kg m −3 was applied as in the case of the Canin East GIP. We determined the mass balance year for balance calculations using a floating-date systems [56], where the surveys were carried out close to the end of the ablation period (see Tables 1 and 2).

Climate Data and Analysis
The meteorological variables (wind fields at 10 m above surface, sea level pressure, precipitation) used in this study in order to analyse synoptic conditions are provided by the ERA-Interim dataset [57]. The ERA-Interim is an atmospheric reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) for the whole globe. The data assimilation system used to produce ERA-Interim includes a 4-dimensional variational analysis (4D-Var) technique with a 12-h analysis window. Additionally, the ERA-Interim dataset is freely available for the period 01.01.1979-31.08.2019 (noting that the ERA-Interim dataset has been superseded by the ERA5 dataset, [58]) at approximately 80 km grid resolution and on 60 vertical levels from surface to 0.1 hPa, which makes this dataset ideal for analysing synoptic situations leading to abundant precipitation or long-lasting rainless events over a specific region in Europe. Figure S1 reports the differences between the aforementioned datasets for winter periods with low-and high snow accumulation.

NAO and AMO Index Data
The North Atlantic Oscillation (NAO) index data and the Atlantic Multidecadal Oscillation (AMO) index data were provided by the National Oceanic and Atmospheric Administration (NOAA). The period between January 1950 and December 2000 for the NAO and 1856 to present for the AMO served as reference for the calculations. Daily and monthly mean NAO index is freely available via NOAA. In this study, we used data based on daily and monthly information. The NAO teleconnection index is obtained through a procedure (Rotated Principal Component Analysis; Barnston and Livezey, 1987), which isolates the primary teleconnection patterns, and is applied to the monthly standardized 500 mb height anomalies over the region embedded between 20 • N-90 • N.

Local Temperature, Precipitation and Glacier Length
Mean summer air temperature (MSAT) and winter precipitation data come from [28] and are updated for the last 6 years from meteorological data recorded at the Canin AWS ( Figure 1) located at an elevation of 2203 m. Snow thickness has been measured since 1972 at the nivological station of Gilberti ( Figure 1) at an altitude of 1830 m a.s.l. (source Regional Administration of Friuli Venezia Giulia, snow and avalanches office http: //www.regione.fvg.it/asp/valanghe/welcome.asp; last accessed 5 August 2020). Data on glacier length come from the Italian Glaciological Committee reports [59]. The sum of melting during the ablation season is computed by using a temperature index, i.e., a degreeday model which is one of the most widely used methods for hydrological and ice-dynamic modelling as well as climate sensitivity studies [47,60]. Although simple, this method is a powerful tool, at times even more reliable than energy balance models at a catchment scale [61]. The amount of melting in mm w.e was calculated by applying a Positive degreeday model with a degree-day factor for melting snow of 4.0 cm day −1 K −1 [47,60,61].

Glacier Area
The area of the Canin East GIP varied substantially over the 12-year observational period (2006-2018) without any decreasing trend ( Table 4). The largest area (27,448 m 2 ) was observed in 2014, which was more than twice (243%) the area of 2006, when the GIP was the smallest in size. The GIP reached the second largest area (19,425 m 2 ) in 2011, with increase of 72% with respect to 2006. The annual areal changes are more or less equivalent around the entire GIP's edge. The GIP commonly constitutes a continuous ice body but has a smaller area at the terminus, which tends to separate from the main ice body during negative mass balance years. The area of the Canin East GIP in 2018 was equal to 15,914 m 2 , which represents 41% increase in area with respect to 2006 (11, Table S1). It is important to note how large fluctuations in ice patch areas are strictly dependent on residual snow and firn of the previous winter(s). Both surveys were performed at the end of the ablation season with no fresh snow on the ground that could have affected the accuracy of data.

Geodetic Mass Balance and Volume Change
The five-year cumulative mass balance of the Canin East GIP was strongly positive with 6. It is worth noting that the calculated annual mass balances of the Canin East GIP have a higher uncertainty than the calculated cumulative mass balances. This is because the GIP surfaces for years 2012, 2014 and 2017 were interpolated using GCPs, whereas the calculations of cumulative mass balances are based on LiDAR LiDAR derived DEMs (Tables 1-3). The uncertainty for annual mass balances can be expressed as the maximum positive and negative elevation difference between LiDAR derived DEM through Delaunay triangulation and GCPs-derived DEM using natural neighbour (Table 3) multiplied by the density of the GIP material, which results in +0.22/−0.28 m w.e.  It is worth noting that the calculated annual mass balances of the Canin East GIP have a higher uncertainty than the calculated cumulative mass balances. This is because the GIP surfaces for years 2012, 2014 and 2017 were interpolated using GCPs, whereas the calculations of cumulative mass balances are based on LiDAR LiDAR derived DEMs (Tables 1-3). The uncertainty for annual mass balances can be expressed as the maximum positive and negative elevation difference between LiDAR derived DEM through Delaunay triangulation and GCPs-derived DEM using natural neighbour (Table 3)   GCP survey, whereas a 5-year and a 12-year surface elevation change reported for the period 2006-2011 and 2006-2018, respectively, refers to Lidar derived surface elevation changes. The GIP outline was adjusted for every surveyed year, and the yearly area change was fully taken into account in the volume change calculations. The large topographic differences observed in some of the surveys, as in the case of the 2012-13, 2013-14 and 2014-15 hydrological years, are due to several different processes able to locally amplify the effect of extreme snow accumulation or ablation: (i) wind snow-blowing; (ii) irregular karst topography due to the presence of karstic shaft, hollows, dolines and canions; (iii) effect of avalanches; (iv) presence of bergschrund.  It is worth noting that the calculated annual mass balances of the Canin East GIP have a higher uncertainty than the calculated cumulative mass balances. This is because the GIP surfaces for years 2012, 2014 and 2017 were interpolated using GCPs, whereas the calculations of cumulative mass balances are based on LiDAR LiDAR derived DEMs (Tables 1-3). The uncertainty for annual mass balances can be expressed as the maximum positive and negative elevation difference between LiDAR derived DEM through Delaunay triangulation and GCPs-derived DEM using natural neighbour (Table 3)  In 2012, when glacier ice outcropped at the end of the ablation season, a set of 7 ablation stakes where drilled by using a Heuke driller in the ice patch ( Figure 2). Positive mass balance in 2013 and 2014 kept ablation stakes covered by snow and firn also during summers 2015, 2016 and until August 2017 when glacier ice outcropped again, and it was therefore possible to measure their topographical position. Only 5 ablation stakes were found, and results show low displacements over the surface of the ice patch with local displacement in the order of 11 to 47 cm in 5 years (i.e., 2 to 10 cm yr −1 ).

Winter Meteorological Patterns in Dry and Wet Winters
During dry winters, the synoptic conditions generally show relatively high sea level pressure (SLP) anomaly, which in fact can be considered as a high-pressure blocking system (or "blocking high"). This is apparent in the cases of 2007-2008 and 2011-12 ( Figure 6). The positive SLP anomaly (reference period: January 1981-January 2010), stable conditions and the downwelling motions together over the southeastern Alps favor the relatively low winter seasonal precipitation. . Mean seasonal precipitation (mm/day, colours) and near-surface (10 m) mean seasonal winds (m/s) are represented along with the mean sea level pressure anomalies (reference period: January 1981-January 2010, unit: hPa). Well-structured and persistent negative sea level pressure anomalies in the Atlantic influenced the regional circulation in Europe and the Mediterranean causing recursive orographic precipitation in the southern Alps, while positive sea level pressure anomalies manifested in low winter precipitation totals. The red square depicts the region of interest.

Climate Data from Local Observations
Analysis of 1979-2018 climatological records (Figure 7) reveals increasing summer (JJA) temperature over the entire period with no signs of a reversal in the trend. Although the interannual variability is unchanged, both annual minima and maxima show a trend towards warmer values. Summer 2003 (11.3 °C) is the warmest observed, followed by summers 2015, 2012 and 2017, which, all together, are the four hottest summers observed since 1851 [28]. The coolest summer of the last 40 years was 1984 at 6.6 °C. Accordingly, the length and magnitude of the ablation period calculated from the PDD model increased. The number of PDD rose from 169 to 187. The potential melting calculated using a degree day factor for melting snow of 4.0 mm day −1 k −1 [60] on average increased about Furthermore, the 2008-2009 winter was characterized by secondary lows forming in the Mediterranean basin with the Alps at the convergence zone of cooler air masses from North and warmer and moist air masses from the Adriatic and the Mediterranean. The combination of the aforementioned conditions created one of the snowiest winters in the Eastern Alps since 1930 [62]. On the contrary, in 2013-2014 winter the subpolar low showed a persistent location of its center slightly south east of Iceland inducing advection of mild and wet south-westerly winds over Europe and the Alps. The long-lasting lowpressure region and a secondary cyclogenesis over the Mediterranean caused recordbreaking winter snow accumulation in the Eastern Alps [63,64]. More specifically, 2% of the winter Vb cyclones exceed the 95th precipitation percentile in the Alpine region [37], which can lead to abundant winter precipitation as occurred in 2013-2014.

Discussion
Although the variability of winter precipitation and glacier accumulation is strongly related to the NAO in Scandinavia and Svalbard, with the AMO also influencing the northernmost glaciers [65,66], this signal over the Alps is somewhat contradictory. Marzeion and Nesje (2012) [67] showed that the NAO in the western Alps is weakly anti-correlated with the mass balance of glaciers due to temperature and winter precipitation anomalies. In the Julian Alps, located in the south-east side of the alpine chain, there is no evidence of correlation between NAO and the past and recent evolution of ice bodies (rs = 0.09; Figure S2). Looking at the long-term evolution of alpine glaciers, Huss et al. (2010) [43] found that glacier mass budget in the Swiss Alps varied in phase with the AMO and is therefore strictly related to North Atlantic variability, but associated with the ongoing global warming. Their results also suggest how the size of glaciers in the Alps reacted in response to AMO-driven variations at least over the last 250 years. They also stated that at the regional alpine scale, AMO is thought to influence the mass budget of glaciers mainly through minima and maxima in temperature of the ablation season.
Regular position measurements of the Canin East glacier terminus from fixed benchmarks started in 1896 [68] and continued through the 20th century up to the present days, with only occasional interruptions during the 1st and 2nd World War [28]. Comparison of glacier-length variation with the AMO index and the normalized MSAT recorded at 2200 m a.s.l. in the Julian Alps over the time window 1851-2018 ( Figure 9) show how multidecadal variations are in phase with AMO with a very strong statistically significance correlation (Sperman's test; rs = −0.94; p < 0.05), whereas the statistic significance with the extended-winter (djfm) 1920-2018 precipitation is weaker (rs = 0.29; p < 0.005). Advances and

Discussion
Although the variability of winter precipitation and glacier accumulation is strongly related to the NAO in Scandinavia and Svalbard, with the AMO also influencing the northernmost glaciers [65,66], this signal over the Alps is somewhat contradictory. Marzeion and Nesje (2012) [67] showed that the NAO in the western Alps is weakly anti-correlated with the mass balance of glaciers due to temperature and winter precipitation anomalies. In the Julian Alps, located in the south-east side of the alpine chain, there is no evidence of correlation between NAO and the past and recent evolution of ice bodies (r s = 0.09; Figure S2). Looking at the long-term evolution of alpine glaciers, Huss et al. (2010) [43] found that glacier mass budget in the Swiss Alps varied in phase with the AMO and is therefore strictly related to North Atlantic variability, but associated with the ongoing global warming. Their results also suggest how the size of glaciers in the Alps reacted in response to AMO-driven variations at least over the last 250 years. They also stated that at the regional alpine scale, AMO is thought to influence the mass budget of glaciers mainly through minima and maxima in temperature of the ablation season.
Regular position measurements of the Canin East glacier terminus from fixed benchmarks started in 1896 [68] and continued through the 20th century up to the present days, with only occasional interruptions during the 1st and 2nd World War [28]. Comparison of glacier-length variation with the AMO index and the normalized MSAT recorded at 2200 m a.s.l. in the Julian Alps over the time window 1851-2018 ( Figure 9) show how multidecadal variations are in phase with AMO with a very strong statistically significance correlation (Sperman's test; r s = −0.94; p < 0.05), whereas the statistic significance with the extended-winter (djfm) 1920-2018 precipitation is weaker (r s = 0.29; p < 0.005). Advances and retreats of the Canin East glacier terminus are in phase with the AMO (Figure 9), which explains almost 120 years of glacier oscillation in the Julian Alps very well. Short periods of glacier advances recorded in 1910-1920 and 1945-1975, as well as in the 1930s, are very well linked and modulated by drops in MSAT. Spearman's correlation test between the two variables returns a moderate to strong correlation (r s = −0.63 p < 0.05). The fast and steady increase in MSAT from early 1980s due to anthropogenic global warming is superimposed to the AMO and explains the dramatic glacier reduction that is particularly evident after the 1980s.
Atmosphere 2021, 12, x FOR PEER REVIEW 17 of 26 retreats of the Canin East glacier terminus are in phase with the AMO (Figure 9), which explains almost 120 years of glacier oscillation in the Julian Alps very well. Short periods of glacier advances recorded in 1910-1920 and 1945-1975, as well as in the 1930s, are very well linked and modulated by drops in MSAT. Spearman's correlation test between the two variables returns a moderate to strong correlation (rs = −0.63 p < 0.05). The fast and steady increase in MSAT from early 1980s due to anthropogenic global warming is superimposed to the AMO and explains the dramatic glacier reduction that is particularly evident after the 1980s. The almost perfect match of glacier reaction to summer temperature forcing, without any apparent delay as in the case of larger alpine valley glaciers, is likely due to the small size of the ice body, which led to a short reaction time (<2 yr) in response to climatic variability [69]. Although both AMO and MSAT explain well the long-term glacier variability in the Julian Alps in agreement with what has been observed in the Alps so far, the recent The almost perfect match of glacier reaction to summer temperature forcing, without any apparent delay as in the case of larger alpine valley glaciers, is likely due to the small size of the ice body, which led to a short reaction time (<2 yr) in response to climatic variability [69]. Although both AMO and MSAT explain well the long-term glacier variability in the Julian Alps in agreement with what has been observed in the Alps so far, the recent drop in envELA and the positive mass balance 2006-2018 are in contrast with what is globally reported elsewhere in the Alps [12].
The general increase in winter (DJFM) precipitation after the negative anomaly of the 1990s (Figure 9a) and the 2008-2009 and 2013-2014 extremes (Figure 8), which brought much higher than average snow accumulation over the ice bodies, are responsible for the observed ice bodies stability and for the 2006-2018 positive mass balance. Recent positive mass balance cannot be attributed to favorable temperature because summers have become longer and hotter (Figure 7d) increasing the snow/ice melting on average by about 29 mm a −1 since 1979 according to the DDM (Figure 7d). Spearman's correlation analysis between the observed 2006-2018 mass balance and the PDD model results shows a moderate correlation statistically not significant (r s = −0.45; p = 0.26). The observed lowering of the envELA in the 2000s, which brought ice body stability and recovery, is therefore due only to changes in precipitation regime and, particularly, in the occurrence of few extreme winter accumulation seasons. This circumstance influenced the fate of ice bodies in the Julian Alps for several years also thanks to the short (<2 yr) metamorphosis of snow into ice. This finding is consistent with Spearman's correlation test between the 2006-2018 annual mass balance and the winter (djfm) precipitation as well as winter (djfma) snowfall, giving a very strong correlation statistically significant (r s = 0.91; p < 0.005 and r s = 0.74; p < 0.05, respectively). Extreme winter accumulation seasons led to multi-annual snow metamorphism and formation of new ice as reported at the end of very hot summers which led to highly negative mass balance in 2012 and 2016. In those years when all previous winter snow was totally removed by melting, bare ice outcropped over the entire surface of the ice patch. The topographic position of ablation stakes indicates a weak dynamics of the Canin East GIP, which is typical of mass movement of ice patches [1] rather than glacier movement due to internal deformation and flow as reported from Montasio glacier by [2,30]. Montasio glacier is the only one to present well defined accumulation and ablation areas (Figure 4), whereas the other ice patches of the Julian Alps do not present a clear ELA. The reason is ascribable to the local topography and the present size of the ice bodies. Snow accumulation due to avalanches affects only the upper part of Montasio glacier, while other Julian Alps ice bodies are often affected by avalanche activity along their entire surface. The dominant role played by solid precipitation in regulating the recent annual mass balance of Montasio glacier as well as the statistically non-significant correlation with temperature of the ablation season was recently also reported by De Marco et al. 2020 [30].
Observations from glaciers in the Austrian Alps (Abermann, Kuhn and Fischer, 2011) confirmed this finding further, suggesting that a slowdown in glacier recession in this area might be related to "above average accumulation". By dividing the Alps into different subregions, each one identified from 2000-2016 annual mass balance time series, Davaze et al. (2020) [15] found that the only climatic parameter statistically significant explaining more than 30% of the observed inter-annual mass balance variance in the Central and Eastern Alps and Tyrol is the precipitable water during winter, that combined with temperature, influences winter precipitation amount.
It has previously been reported how the maritime glaciers in the Julian Alps recently showed higher climatic sensitivity to precipitation than to air temperature (Colucci and Guglielmin, 2015) (Salvatore et al., 2015). Other authors recently highlighted how residual ice masses in other maritime areas (i.e., where MAP is high) in the southern Alps showed contradictory signs of stability and/or shrinking slow down instead of rapidly disappearing as might be expected (Scotti and Brardinoni, 2018).
Topoclimatic factors and the role of wind-blown and avalanche-fed accumulation (Fischer, 2011 Explaining the long-term decadal variability of winter precipitation in the Julian Alps with AMO is not an easy task, both because the sample is not statistically significant (r s = −0.12; p = 0.26) and the signal not as clear as for summer temperature, as already reported by Napoli et al. (2019) [44]. Nevertheless, several authors recently stated that the intensity of heavy precipitation events seems to increase more than mean precipitation under a warmer climate in most of the Euro-Mediterranean region [70,71] consistently with a larger moisture-holding capacity and the Clausius-Clapeyron law. Giorgi et al. (2011Giorgi et al. ( . 2016 [72,73] highlight the role of topography in considerably modulating the precipitation change signals in global climate projections, with orographic convection and shading being the most important factors. Their projections for the summer season show an increase in convective rainfall due to enhanced potential instability by high-elevation surface heating and moistening. They also stated how this should also affect the occurrence of extreme events. The importance of orography in modulating precipitation variability in the Alps is highlighted by Napoli et al. (2019) [44]. They noticed how the distribution of precipitation has recently varied over time, with generally higher precipitation at high altitude compared to low elevations, although a straightforward definition of the main factors driving such behavior is not possible at present knowledge.
The three elements involved in orographic precipitation are (i) moist and large-scale flow towards an obstacle; (ii) mesoscale orographic lifting (which cools the air to saturation) and condensation; (iii) conversion of condensate to precipitable particles involving microphysics of the clouds, smaller scale convection and turbulence, which in turn locally enhance the growth and fallout of precipitation [74]. Increase in precipitation is a sensitive function of low-level thermodynamic structure (the planetary boundary layer, PBL), and processes may differ depending on whether or not the low level flow undergoes blocking (flow around) or is unblocked (flow over) [75]. Parameters controlling this two outcomes are the wind speed impacting the mountain chain, the height of the chain and the degree of stability of the air mass measured by the Brunt-Väisälä frequency [74]. In the case of flow over, a higher liquid content in the clouds generated by the orographic-forced updraught allows liquid water to attach to precipitation particles by mechanisms of coalescence and riming, further growing the snow flakes and producing the most intense precipitation over the first reliefs of the chain [76].
Recently, several studies also focused on the impact on heavy precipitation caused by interaction between the sea and the atmosphere (e.g., [77]) pointing out the crucial role of sea surface temperature (SST) in influencing turbulent heat and water vapor surface fluxes (e.g., [78,79]). Particularly, the moist southeasterly wind (Sirocco) from the Adriatic is often responsible for long lasting orographic precipitation in the Pre-Alpine and Alpine areas of NE Italy. In case of conditionally unstable flow, the interaction with the orography eventually led to the formation of intense convective system upstream of the orography in all seasons [80][81][82]. Stocchi and Davolio [83,84] have shown how during Sirocco events, SST plays an indirect role in determining the location of precipitation, influencing the characteristics of the PBL and the low-level jet as well as its interaction with the orography. The atmospheric water budget is not noticeably influenced by sea surface evaporation regardless of higher SST, which is instead primarily modulated by synoptic mesoscale patterns. On the other hand, higher SST influences the thermodynamic structure of the PBL determining a more likely transition from flow around to flow over conditions. The consequence of increased SST is to intensify vertical motions and to shift towards the mountain crest downstream the area with intense orographic uplift and condensation.
It is worth noting that SST in the Adriatic Sea are lately increasing at accelerated pace, i.e., from +1.3 • C per century in the 1946-2015 period [85] to +3.76 • C per century in the 1979-2015 period, passing from a late 1970s winter average of 8.89 • C to 10.28 • C in 2015 [85]. SST in the North Adriatic is subjected to high variability amplified by a shallow basin due to changes in wind direction associated with synoptic conditions. Such variability should be detailed and analyzed on a daily basis in association with each precipitation event, but at this stage this goes outside the scope of this paper. Nevertheless, Spearman's correlation rank analysis between cumulated winter (djf) precipitation ( Figure 6) and mean winter (djf) SST returns a statistically significant direct association between the two variables over the period 1979-2015 (r s = 0.325; p < 0.05).
It is likely that for all the aforementioned evidence Adriatic SST also contributed to the exceptional winter snowfalls recorded in the 2000s. Analysis of wet winters presented in this work showed how weather patterns where favorable to recurrent conditions able to trigger Scirocco events and thus orographic precipitation. In fact, both in the winter 2008-2009 and 2013-2014, a well-structured and persistent negative sea level pressure anomalies influenced the regional circulation in Europe and the Mediterranean (Figure 10).
Atmosphere 2021, 12, x FOR PEER REVIEW 20 of 26 [85]. SST in the North Adriatic is subjected to high variability amplified by a shallow basin due to changes in wind direction associated with synoptic conditions. Such variability should be detailed and analyzed on a daily basis in association with each precipitation event, but at this stage this goes outside the scope of this paper. Nevertheless, Spearman's correlation rank analysis between cumulated winter (djf) precipitation ( Figure 6) and mean winter (djf) SST returns a statistically significant direct association between the two variables over the period 1979-2015 (rs = 0.325; p < 0.05). It is likely that for all the aforementioned evidence Adriatic SST also contributed to the exceptional winter snowfalls recorded in the 2000s. Analysis of wet winters presented in this work showed how weather patterns where favorable to recurrent conditions able to trigger Scirocco events and thus orographic precipitation. In fact, both in the winter 2008-2009 and 2013-2014, a well-structured and persistent negative sea level pressure anomalies influenced the regional circulation in Europe and the Mediterranean (Figure 10). In 2008-2009 recursive lows formed in the Mediterranean with the Alps at the convergence of colder northerly fluxes and the milder and wetter south flows, while in 2013-2014 a strong negative anomaly of the semi-permanent sub-polar low repeatedly triggered warm and wet winds towards the southern side of the Alps producing heavy precipitation. The location and duration of such events at mid-latitudes are largely controlled by atmospheric circulation where the jet-stream strongly meanders and the eastward transport of weather systems slows down [86]. Meanders in the jet are referred to as Rossby waves, recognizable from the meridional wind anomalies at 300 hPa [87], and when they are slower and more amplified, the occurrence of extreme weather conditions in terms of high anomalies of land surface temperature and land precipitation is favored [88].
During episodes of amplified Rossby waves, circulation patterns at upper and lower pressure levels are vertically aligned and precipitation extremes generally occur below a trough, thus in cyclonic conditions [87]. In the winters of both 2008-2009 and 2013-2014, mid-latitudinal upper tropospheric circulation was characterized by a strongly meandering jet which especially in 2013-2014 encircled the Northern Hemisphere in a regular wavenumber 4 pattern (Figure 10f,n). This configuration gave conditions for recursive and frequent precipitation in the Southern Alps opposite to conditions observed during seasons 2007-08 and 2011-12 when meridional wind anomaly was strongly negative favoring dry conditions in the southern Alps (Figure 10c,i). In both cases, the source of moisture reaching the southern Alps below the PBL was the Adriatic Sea.
The occurrence of particularly snowy winters in 2008-2009 and 2013-2014 providing accumulation for perennial ice and firn patches that were able to survive very warm summers is also reported by Scotti and Brardinoni (2018) [29] from two small cirque glaciers in the Orobie Alps, central Southern Prealps. This area shares several common features with the Julian Alps in terms of orographic precipitation, and common aforementioned processes are involved in producing flow around conditions during Scirocco events, including the main moisture contribution in the PBL from the Adriatic Sea [44,74]. Different location of the center of the low pressure systems resulted in generally close-to-or lower-than average winter temperatures in the Alps in 2008-2009 in contrast with the strong positive temperature anomaly of 2013-2014. On the other hand, positive temperature anomalies over the Arctic were recorded in both winters (Figure 10b,e) as well as positive sea level pressure anomaly, more pronounced in 2008-2009 over Alaska, Siberia and the North Pacific range. The past few decades have been characterized by Arctic warming, which has almost doubled compared with the entire northern hemisphere due to a combination of increased anthropogenic global warming and positive feedbacks leading to the so-called Arctic Amplification (AA) [89,90]. Observations and modeling recognized several large-scale modifications in the atmospheric circulation in association with sea-ice loss and earlier snow melt, which in turn affect precipitation, seasonal temperatures, storm tracks and surface winds in mid-latitudes [91]. Singular strong weather events normally have a dynamical origin and many of them result from persistent weather patterns, typically associated with blocking and high amplitude waves in the upper-level flow. Evidence links AA with an increased tendency for a slower eastward movement of Rossby Waves that in turn favor weather extremes triggered by persistent weather conditions due to amplified flow trajectories [91]. The persistent weather conditions associated with severe snowfalls in the southern Alps in 2008-2009 and 2013-2014 exacerbated by local orography and increasing SST, which led to positive glacier mass balance in few sectors of the Alps, remind of the Brikdalsbreen event on a smaller scale. This event is a winter precipitation-induced decadal-scale glacial advance observed in southern Norway in the 1990s due to recurrent positive NAO [65].

Conclusions
The 12-year positive mass balance of small and residual ice bodies (very small glaciers and ice patches of glacial origin) in the Julian Alps, southeastern European Alps, in the period 2006-2018, provides the opportunity to analyze the contingent weather patterns driv-ing such anomalous situations. While multidecadal variability of AMO, which modulates the temperature of the ablation season, explains 120 years of glacier advances and retreats, recent stability and recovering of minor glacial bodies in few sectors of the Alps as well as lowering of the envELA need different explanations. Areas characterized by the occurrence of intense orographic precipitation under certain synoptic circumstances, as in the case of few sectors of the southern Alps, recently recorded extreme and prolonged winter snowfalls that are sufficient to offset longer and warmer summers. Such extreme seasonal events are favored by the thermodynamic structure in the boundary layer, likely accentuated by accelerated increase in the Adriatic SST. Several extreme weather events, which generally have a dynamical origin, are initiated by persistent weather patterns associated with slower eastward movement of planetary waves (Rossby waves) having more meandering paths. Evidence links the occurrence of amplified upper-level flow trajectories to the Arctic Amplification and, therefore, to ongoing global warming. Although further summer warming is expected in the next decades, modification of weather patterns due to a set of teleconnections and higher occurrence of extreme and persistent snowfall events in winter might represent a crucial input in the survival of small glacial remnants in maritime mountain areas in the near future.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4 433/12/2/263/s1, Figure S1. Seasonal precipitation based on the ERA-Interim (left column) and ERA5 datasets (right column) during winters (December-January-February) with relatively low snow accumulation (2007/2008 and 2011/2012 first row and second row, respectively). The red square depicts the region of interest. Units are mm/day. Figure S2. NAO index vs winter precipitation divided in winter Snow Total Accumulation (STA) and Winter (djfm) Precipitation (WP). STA is the sum of all single snowfalls manually measured on the ground between the second half of November and the end of April. Association between WP and NAO (rs = −0.30078, p = 0.05931; R2 = 0.1754) is statistically not significant. Table S1. Area, volume change and mass balance of all ice bodies in the Italian Julian Alps for the period 2006-2018.
Author Contributions: R.R.C. initiated this research, led the team and wrote the manuscript together with collecting contributions from the other authors. M.Ž. analysed the annual and long-term geodetic mass balance and prepared related maps as well as wrote the corresponding chapters. C.D.G. set the degree day model calculation and took part in most of the surveys on Canin GIP. E.M. elaborated the LiDAR information and registered the digital elevation models writing the corresponding chapter. S.P. collected the GPS glacier monitoring and elaborated the raw data. C.Z.T. worked on climate data and analysed the weather patterns. N.F.G. helped frame the sections on glacier response to climate and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.