The glaciers of the Dolomites: the last 40 years of melting

Small Alpine glaciers located below the regional equilibrium line altitude are experiencing considerable ice loss and are expected to fragment into smaller glacial bodies and eventually disappear. Monitoring such glaciers by satellite remote sensing is often challenging because their size and surrounding topography are incompatible with the current spatial resolution of non-commercial satellites. The Italian Dolomites (southeastern Alps) are a region clearly illustrating such challenges and where no long-term glacier mass balance data are available. This renowned Alpine sector hosted tens of glaciers up until a few decades ago, with now only 12 remaining. This study presents a multi-decadal (1980s–2023) estimation of surface elevation change and geodetic mass balance of the current mountain glaciers present in the area. Calculations are based on geodetic data: high resolution and accuracy are obtained with uncrewed aerial vehicle (UAV), structure from motion (SfM), and airborne light detection and ranging (lidar) methods, from 2010 to 2023. SfM on historical aerial imagery is used for previous decades. We found an average cumulative surface elevation change of $-\mathrm{28.7}\pm \mathrm{2.6}$ m from 1980s to 2023, 33 % of which occurred between 2010–2023. The average geodetic mass balance rate for the whole period is $-\mathrm{0.64}\pm \mathrm{0.06}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (water equivalent), varies widely among sites, and is less negative than the Alpine reference glaciers' mass balance. Regionally, 66 % of the volume loss is related to the Marmolada glacier alone. Losses of ice mass and area show that the Dolomites are rapidly losing their glaciers.

1 Introduction

Glaciers worldwide have been losing mass at alarming rates over the past decades (Zemp et al., 2019; Hugonnet et al., 2021). This is particularly evident in regions where warming is occurring at a faster rate than the global average (Rantanen et al., 2022; ICCI, 2022). Among these regions are the European Alps (Hock et al., 2019), where severe impacts on the cryosphere have been observed because of rising temperature (Auer et al., 2007), changes in seasonal precipitation patterns, total radiation, humidity, and snow accumulation (Gobiet et al., 2014; Rumpf et al., 2022).

The European Alps can effectively be considered the birthplace of modern glaciology (Clarke, 1987). Despite the progressive improvement of in situ (Huss et al., 2015; WGMS, 2021) and remotely sensed cryosphere monitoring efforts for the European Alps (Davaze et al., 2020; Paul et al., 2020; Sommer et al., 2020), areas remain where only scattered data and observations are available. One of those regions is the Italian Dolomites (southeastern Alps).

The Dolomites are mostly known for their landscapes, resulting from an extremely rich geological diversity (Panizza, 2009). They represent a complex area from a socio-economic, geographical, and cultural point of view, and for this reason they have been studied intensively in terms of geology, geomorphology, and biodiversity (Bosellini et al., 2003; Panizza, 2009; Pignatti and Pignatti, 2014). However, glaciers of this region have always remained on the margins of scientific exploration.

Glaciers in the Dolomite region were numerous until a few decades ago (CGI-CNR, 1962). The complex topography of the area ensures that many sites are orographically protected from direct sunlight and favour the accumulation of snow through avalanches. Until recently, such a context allowed for the existence of glaciers well below the environmental equilibrium line altitude (envELA) (Castiglioni, 1925; Žebre et al., 2021), which in the Eastern Alps during the 20th century ranged from 2700 to 3100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (above sea level) (Žebre et al., 2021). Only a few massifs had an accumulation basin located above the envELA, allowing for larger ice bodies to develop, such as the Marmolada and Fradusta glaciers.

Assessing the state of glaciers through their mass balance is considered an essential component in the global-climate-related set of observations (Zemp et al., 2009). At the regional scale, this involves generating large composite datasets through in situ measurements and remote sensing campaigns. Various techniques have been developed over the last few decades to improve remote measurements of glacier mass change (Berthier et al., 2023). Workflows involving sensors mounted on uncrewed aerial vehicles (UAVs) or helicopters are particularly effective for monitoring smaller topography-bounded glaciers. Among such settings, one technique that has been successfully used is structure from motion (SfM) (Carrivick et al., 2016). SfM allows for both terrestrial (e.g. Piermattei et al., 2015; Marcer et al., 2017) and aerial (e.g. Smith et al., 2016) high-resolution surveys relying on photos taken from different positions. The same workflow can be applied to historical imagery initially acquired for mapping purposes, allowing for the retrieval of past glacier change (Mertes et al., 2017; Knuth et al., 2023). Despite the large improvements in the applications of SfM, the state of the art for high-resolution topographic surveys remains light detection and ranging (lidar) (Bhardwaj et al., 2016). Lidar is an active remote sensing technique that uses a laser pulse and its two-way travel time to reconstruct the surrounding scene. It can be used to monitor surface elevation changes at different scales, depending on the sensor and platform used (e.g. Knoll and Kerschner, 2009; Fischer et al., 2016; Okyay et al., 2019; Securo et al., 2022).

To the best of our knowledge, no specific studies have been conducted regarding the current status of Dolomite glaciers or their decadal change. The aim of our work is to fill this gap, providing a description of the mountain glaciers in the region and quantifying their evolution during the last 40 years. To achieve these objectives, surface elevation changes and the geodetic mass balance of the 12 current Dolomite mountain glaciers were calculated from the 1980s to 2023, with almost one measurement per decade. Available temperature, precipitation, and snow accumulation data from the regional network of weather stations are used to contextualise the results.

2 Previous glaciological research in the Dolomites

The first study to mention the Dolomite glaciers dates to 1888 (Richter, 1888). In his volume, the Austrian glaciologist described the then-known glaciers in the eastern Alpine region, including two of the largest ice bodies of the Dolomites: the Marmolada and Travignolo glaciers. However, the first work introducing topographical maps and quantitative information was compiled by Olinto Marinelli (Marinelli, 1910). Thereafter, the available information about the glaciers of the Dolomites mostly comes from works that have dealt generally with glaciers in the Italian Alps (e.g. Porro, 1925; CGI-CNR, 1962; Smiraglia et al., 2015).

Data presented for the wider Italian Alps or pan-Alpine products are not always accurate for the Dolomites. In the past, only sporadic attention, from a glaciologic point of view, has been paid to specific massifs (Castiglioni, 1925, 1930; Nangeroni, 1938; Del Longo et al., 2001; Cibien et al., 2007). Some more information is available from institutional reports edited for administrative regions, but considering the incompleteness of these dispatches, they are not discussed here. According to the Italian Glacier Inventory (CGI-CNR, 1962), in the late 1950s, the Dolomites hosted 33 glaciers, with a total surface of 8.19 km².

The last partial survey of the Dolomite glacier area was carried out by the Agency for Environmental Prevention and Protection of Veneto (ARPAV), which analysed the surface variations for the 27 main glacial bodies of the Dolomites (Crepaz et al., 2013), representing the 72 % of the total glacierised area in the region. Results show an area loss of approximately 50 % from 1910 to 2009. The most recent inventory available for Italian glaciers (Smiraglia et al., 2015) reports that 51 glacial bodies were present in the Dolomites in 2009, spanning a total area of 5.04 km² and making up 1.4 % of the total Italian glacierised area. Among the 51 glacial bodies, 13 are classified as mountain glaciers, while 38 are considered snow or ice patches. When we use the term ice patch, we refer to the description of ice patch of glacial origin proposed in Serrano et al. (2011).

The only glacier of the Dolomites that has been extensively studied in the past years is the largest in the area, the Marmolada Principale glacier. It was surveyed with ground-penetrating radar in 2004 and 2014, revealing a 30 % reduction in thickness, hinting at an imminent fragmentation and a likely complete disappearance before 2050 (Santin et al., 2019). The transition between glacial and periglacial landforms in the area is discussed in Seppi et al. (2014), while studies focused on rock glaciers (Krainer et al., 2010, 2012) and debris-covered glaciers (Securo et al., 2024a) have also been done recently. Despite their limited size, the glaciers and glacial bodies in the Dolomites are also investigated in relation to geohazards. This is firstly confirmed by the deadly Marmolada ice and rock avalanche that occurred on 3 July 2022 (Olivieri and Bettanini, 2023; Bondesan and Francese, 2023) but also by debris flow events on Monte Pelmo, closely related to the Val d'Arcia glacial body (Del Longo et al., 2001; Chiarle et al., 2007), and by the instability of moraines from the Little Ice Age (Zanoner et al., 2017).

As this work specifically focuses on the remaining mountain glaciers of the Dolomite area, attention is given to the following ranges (Fig. 1): Popera, Cristallo, Sorapiss, Antelao, Marmolada, and Pale di San Martino. Other Dolomite massifs still host minor ice deposits, but now being devoid of any dynamics (Smiraglia et al., 2015), they are not considered in this work. The Fradusta and Marmolada glaciers, which recently experienced fragmentation, will still be treated as single glaciers (Table S1 in the Supplement).

Figure 1 (a) The Dolomites in the European Alps (© Map Tiler; © OpenStreetMap contributors 2023, distributed under the Open Data Commons Open Database License (ODbL) v1.0). (b) Location of the six mountain areas hosting active glaciers in the region; location of snow ice patches (white dots), mountain glaciers (light-blue dots) (Smiraglia et al., 2015), automatic weather stations (AWSs) (yellow), and snow depth observations stations (red). (c) 3D views of the Dolomite glaciers derived from 2014 lidar data and IDs of the glaciers from the Italian Glacier Inventory (CGI) (Smiraglia et al., 2015). A few AWSs used in this study are located outside of the map in (b).

3 Data and methods

In this study aerial photographs spanning 3 decades (1980–2012), airborne lidar (2010–2014), and UAV and aerial surveys (2023) were used to produce dense point clouds of the Dolomite glaciers and assess their change through time. Data acquisition and pre-processing, point cloud generation, point cloud alignment and error assessment, and weather station data processing are described in four separate sections (Fig. S1 in the Supplement). When decades such as the 1980s, 1990s, and 2000s are mentioned in the text, we are always referring to multiple years of available data (Fig. 1).

3.1 Acquisition and pre-processing

All imagery used was pre-processed using Adobe Photoshop (v2023) to improve the exposure of heavily shadowed areas. The archive photos span from the 1980s to 2012 (Table S2 and Fig. S2). Digital photos from 2010 and 2012 have an optimal quality for SfM processing but lack metadata. Older analogue aerial imagery was scanned with a non-photogrammetric scanner and therefore lacks distortion and internal orientation information that should be used for camera calibration.

Airborne lidar surveys were performed in 2010 and 2014 using an Optech ALTM 3100 instrument mounted on a helicopter. The produced dense point clouds included areas outside of the glaciers. The ground control points (GCPs) used during georeferencing of all aerial photos were retrieved manually from 2010 lidar point clouds using the tool for point list picking in the open-source software CloudCompare (v2.12). Prior to GCP retrieval, the PCV (portion de ciel visible) plug-in, a generalisation of Tarini et al. (2003) algorithms, was used to improve 3D visualisation. Preference was given in picking points that were steady and well recognisable over time also in older photos, such as big boulders and bedrock features. The number of GCPs per area varied between 6 and 30 depending on glacier size and image quality. The 2023 UAV and aerial surveys lacked GCPs collected on site and therefore needed the same use of lidar-derived points to improve accuracy and precision.

3.2 Structure from motion (SfM) workflow and area mapping

This work used structure from motion (SfM) to retrieve the former glacier surface as dense point clouds, which were later compared with 2010 lidar data, taken as the most reliable reference. The SfM workflow was entirely performed using the software Agisoft Metashape (v2.0.4) and is similar to part of the processing of photos by Knuth et al. (2023) that have no camera calibration information or fiducial marks available. The archive images were adjusted and standardised by manually picking pseudo-fiducial markers (i.e. a consistent point in the photo frame notches) from each photo frame and performing the calibration in Agisoft Metashape. The subsequent photo alignment was done using 200 000 tie points and 20 000 key points with the maximum quality, generating for each scene a sparse point cloud. GCPs were manually placed and adjusted in all the photos.

GCP known-coordinate input allowed for the adjustment of the sparse-point-cloud position, scale, camera positions, and subsequent dense-point-cloud generation. A few scenes also required the computation of a triangular mesh generated by linear interpolation of points to close data voids in the dense point clouds. The same mesh with additional texture can be used for visual interpretation of results (Fig. 2). Digital elevation models (DEMs) and orthophotos were also generated for each survey to enable the manual mapping of glacier areas during all reconstructed time steps, done in QGIS (v3.28). We used the surface lowering observed during the multi-temporal comparisons combined with direct observations on site to assess the glaciers' end when the presence of debris cover prevented the delineation of the glaciers' perimeter using mapping alone. The report of SfM output is available in Table S3.

Figure 2 Examples of the reconstructed 3D textured meshes of the Antelao Superiore, Fradusta, and Marmolada Principale glaciers during different decades. Reconstructions of the Popera, Sorapiss Occidentale, Cristallo, Antelao Inferiore, and Travignolo glaciers are shown in Fig. S3.

3.3 Point cloud comparison and error assessment

SfM and lidar dense point clouds were used to detect and quantify surface elevation change. All of the processing was conducted in CloudCompare (v2.12). The 2010 lidar point clouds are the most reliable product available in this study and were used to refine alignment through manual point picking, adjusting the misaligned 1980s, 1990s, and 2000s point clouds. Once aligned, point clouds from different years were compared to 2010 lidar using the Multiscale Model to Model Cloud Comparison (M3C2; Lague et al., 2013) algorithm. M3C2 is a robust solution to assess the differences between point clouds without the need for surface meshing or DEM generation. The M3C2 workflow efficiently manages complex terrain features like flat and vertical surfaces within a single scene and was already used in similar applications (e.g. Midgley and Tonkin, 2017; Mishra et al., 2022). In certain reconstructions of the older scenes (1980s), the point cloud output of solely SfM contained small void areas in correspondence with the underexposed sectors or because of snow. In those cases we used the Metashape-produced 3D triangular meshes that were sampled and processed in CloudCompare as the other point clouds thereafter.

The M3C2 output that holds the computed differences was rasterised and imported in QGIS (v3.28). Here, we used the mapped polygons of the glaciers during all available time steps to compute zonal statistics. The statistics were computed on two different polygons for each site: common area (i.e. the smaller and more recent outline of the glacier, in times of glacier retreat) and total area (i.e. largest glacier area between the two time steps of each comparison). A density conversion factor of 850 kg m⁻³ (Huss, 2013) was used to convert surface elevation change values of the common area to geodetic mass balance. Each comparison included 2010 lidar data as a reference, and the comparisons were never done using two historical SfM point clouds at a time to reduce sources of uncertainties.

The error in surface elevation change varies for each comparison and comes from different sources: (i) lidar vertical accuracy (δ_L), (ii) alignment error between the point clouds on stable terrain (δ_AL), and (iii) measurement error (i.e. distance uncertainty, δ_M3C2). δ_L was quantified by the lidar surveyors that compared the dataset's vertical accuracy with known points measured with differential GPS and was ±0.12 m. δ_AL was measured for each individual comparison using M3C2 in the frontal and lateral areas outside of the glaciers that share an average slope value to the latter. Slope similarity between glaciers and the stable terrain was calculated in CloudCompare using the geometric feature verticality (Harshit et al., 2022), computed on a cylinder with a search radius of 5 m for all point clouds. Areas of the point clouds above and below the mean glacier verticality ±2σ (95 % confidence interval) and inside the glaciers area were filtered out before using M3C2 again, this time for δ_AL calculation. δ_M3C2 was available as a direct output of the previous processing (i.e. distance uncertainty), and its mean value per glacier was retrieved using zonal statistics. δ_L, δ_AL, and δ_M3C2 were summed in quadrature to estimate the overall uncertainty in surface elevation change (δ_S). The contribution of uncertainties in the area was not considered during this phase.

The error in geodetic mass balance (δ_MB) estimations included an additional uncertainty interval in density conversion (δ_ρ) of $\pm \mathrm{60}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (Huss, 2013). This was added to the δ_S using the product error propagation formula of Eq. (1) (Taylor, 1997):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{\mathit{\delta }}_{\text{MB}}=\left|S\mathit{\rho }\right|\sqrt{{\left({\displaystyle \frac{{\mathit{\delta }}_S}{S}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\delta }}_{\mathit{\rho }}}{\mathit{\rho }}}\right)}^{\mathrm{2}}}.\end{array}$$

Considering our dataset, the reconstruction of older glaciers' surface using scanned imagery was limited in accuracy because of (i) a lack of photogrammetric quality scans of the negatives, (ii) the presence of heavily shaded zones or clean snow, and (iii) the reduced availability of suitable spots for GCP picking in some scenes. These sources of errors are harder to quantify and reduced the number of historical photos suitable for such an application. Further sources of uncertainty may be related to the temporal coverage of aerial photography that was not uniform over all glacial bodies and years (Table S2).

3.4 Weather station network

Precipitation and temperature time series from 20 automatic weather stations (AWSs) of the ARPAV network were integrated with the datasets from two AWSs managed by the autonomous province of Bolzano, which altogether guarantee an adequate spatial coverage of the study area (Fig. 1). The selection of the AWSs was based both on their spatial distribution and proximity to the six mountain areas hosting Dolomite glaciers, as well as on the availability of continuous time series of temperature and precipitation starting at least 20 years before 2022. A linear interpolation was employed to fill temperature gaps of 3 d or fewer, while no reconstruction was performed on precipitation data. Additionally, years with missing data exceeding 5 % of the accumulation (November to April) or ablation (June to August) season were excluded from the analysis. This was implemented at the level of individual AWSs, ensuring the availability of data for each year after averaging across all stations. All the time series begin between 1985 and 2001 and end between 2020 and 2022, with 77 % of them spanning over more than 30 years and 9 % spanning less than 26 years. Due to fewer AWSs being deployed in the 1980s and as a consequence of data cleaning, only 36 % to 64 % of the stations present temperature (23 %–64 % for precipitation) data from 1985 to 1991. However, from 1992 onwards, the data coverage increases drastically.

For each station the standardised anomaly index (SAI) (Katz and Glantz, 1986) was calculated, defined in Eq. (2):

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\text{SAI}={\displaystyle \frac{x_{\mathrm{a}}-x_{\mathrm{m}}}{\mathit{\sigma }}},\end{array}$$

where x_a is either the total precipitation during the accumulation season (for the precipitation SAI, Pr SAI) or the mean temperature during the ablation season (for the temperature SAI, T SAI). x_m and σ are the mean and standard deviation of temperature and precipitation over the reference period 1991–2020. The accumulation and ablation seasons were defined according to the local climatology, taking into account that on average snow cover increases between November and April and decreases from May onward. Finally, SAI values were averaged across all weather stations, resulting in unique Pr and T SAI time series representing the entire region. The pre-processing applied to AWS data may result in an underestimation of total precipitation and therefore of the Pr SAI. Conversely, with regards to temperatures, we expect that the reconstructed data do not significantly alter the mean, as they are prone to distribute randomly both above and below the average.

Daily snow depth measurements from four stations managed by ARPAV provided snowfall information between the 1980s–1990s and 2023, occurring above 1900 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1). These measurements were obtained from both automatic and manual snow monitoring stations. The automatic stations deploy an ultrasound sensor and record a value every 30 min, while the manual measurements are taken once per day at 08:00. The Fedaia station has the longest time series, starting in 1980, while Ra Vales, Monti Alti di Ornella, and Col dei Baldi began recording data in 1993, 1986, and 1987, respectively. Using these data, we show the October–June snow depth on the ground for the most relevant years of our study (1982, 1992, 2010, 2014, 2023). Additionally, we calculated the October–June snow depth on the ground averaged over the five reference years for each station, as well as the total annual snow accumulation from the 1980s to 2023.

4 Results

The dataset used in this study allowed for obtaining the topographic surface and area changes in Dolomite glaciers from the 1980s to 2023, enabling the estimation of their geodetic mass balance. Regarding the nine glaciers analysed, the surface of Sorapiss Occidentale has been reconstructed since 1980; that of Antelao Superiore, Antelao Inferiore, Marmolada, Travignolo, and Fradusta has been reconstructed since 1982; and that of Popera Alto, Popera Pensile, and Cristallo has been reconstructed since 1992 (Figs. 2 and S3).

4.1 Area loss

The manually mapped areas of the Dolomite glaciers (Table S4) showed a retreat from 4.11 to 1.81 km² between the 1980s and 2023 (−56 %). The loss during the last 13 years (2010–2023) was 0.88 km² (Fig. 3), corresponding to a relative reduction of 33 %. The area losses presented here only refer to the glaciers selected in this study and as such are an underestimation of the total area loss involving all the ice bodies in the Dolomites. In the 1980s and 1990s the Dolomite glaciers were more numerous, with several of them having completely melted, having turned into permanent ice patches without apparent ice dynamics, or having been buried by debris.

Figure 3 Areas of the Dolomite glaciers in 1980–1982 (dashed line), 2010 (grey), and 2023 (blue) and debris cover in 2023 (dots). The oldest year available for Popera and Cristallo is 1992 (*). Cristallo 2023 debris cover data are not available due to the presence of snow of negligible thickness at the time of acquisition.

Relative area loss (1980s–2023) showed a wide range of values, between 9 % and 89 %. Where surface debris cover was abundant like in Popera Alto, Travignolo, and Sorapiss, we observed only minor area losses, especially during the last decade. The glaciers less surrounded by steep topography were the two affected by the largest area change from the 1980s: Fradusta (−89 %) and Marmolada (−60 %). The latter were also the only glaciers that underwent fragmentation, with the lower part of Fradusta (F. Inferiore) being the only ice body completely extinct among our dataset. All area variations were almost completely related to glacier fronts except for Marmolada, where even the upper portions were affected by retreat.

The Marmolada glacier still remains the largest glacier of the Dolomites, contributing to 66 % of total area losses from 1980s and currently representing 55 % of the total glacier area of the Dolomites.

4.2 Surface elevation change and mass balance

The average surface elevation change per site was calculated for three periods: 1980s with −5.2 m, 1990s–2010s with −14.1 m, and 2010s with −9.3 m. Due to the lack of data in 1999 and 2001 it was not possible to resolve the 1990s–2000s period at decadal frequency (Table S5).

Between 1980s and 2010 the Dolomite glaciers experienced the greatest reduction in thickness in their frontal parts, except for Popera Pensile and Sorapiss Occidentale, which instead showed greater losses in their upper portions (Fig. 4). The Marmolada, Fradusta, and Antelao Superiore glaciers showed the highest cumulative thickness losses, exceeding 50 m at their former fronts.

Figure 4 Surface elevation change (m) of the Dolomite glaciers from 1982 (1980 for Sorapiss, 1992 for Popera and Cristallo) to 2010. Average surface elevation change is reported for every glacier using the common area (red) and total area (brown). Surface elevation change uncertainty (δ_S) is reported for each comparison in grey.

The largest values of glacier surface elevation loss between 1980s and 2023 can be observed in the Fradusta glacier, with an average loss of 49.70 m. The smallest changes occurred in the Popera Pensile and Popera Alto glaciers, with −17.0 and −19.9 m, respectively (Table S5).

Surface elevation change calculations during the 2010–2023 interval were more accurate than for the 1980s–2010 period, with an uncertainty δ_S ranging between 0.3–0.4 m, and allowed for a uniform comparison of the glacier behaviour during the last decade. Overall, δ_S measurements do not show an evident relationship with slope and have a high variability between sites (Fig. S4), influenced by the uneven distribution of samples across slope bins. During 2010–2023 the acquisition was consistent with dates close to the end of the ablation season, and all glaciers were surveyed (Fig. 5a). The average surface elevation changes between 2010 and 2014 are shown in Fig. S5, with five glaciers (Popera Alto and Pensile, Sorapiss, Antelao Inferiore, and Marmolada) that recorded a surface elevation increase since 2010. During 2010–2023 Marmolada Principale, Marmolada Ovest, and Antelao Superiore experienced greater losses, exceeding 30 m. The highest absolute values, corresponding to almost −35 m, were reached in the area involved in the ice avalanche that happened in a detached part of Marmolada Principale, on 3 July 2022 (Olivieri and Bettanini, 2023; Bondesan and Francese, 2023). These are shown by the kernel density plots of surface elevation loss (Fig. 5b, Marmolada Collapse). The Fradusta Inferiore glacier was not included in the common area measurements as it had already disappeared before 2023 surveys took place.

Figure 5 Surface elevation change (m) measured across the Dolomite glaciers from 2010 to 2023. (a) Average surface elevation change is reported for every glacier using common area and total area. Surface elevation change uncertainty (δ_S) is reported for each comparison in grey. (b) Kernel density distribution plot for surface elevation loss with mean loss (dotted black line), with all former Marmolada glacier sections presented as divided plots.

During the last 13 years we observed the only thickness increases in the central part of the Sorapiss Occidentale glacier, which is currently debris covered (Fig. 3). On that glacier a rise of more than 10 m was observed close to a wide serac whose presence is possibly related to a small surge induced by a recent rockfall (Fig. 5a) in the accumulation area as well as by internal glacier dynamics. Rockfall debris and boulders may also have contributed to the increase in surface elevation. A very small area of the Marmolada Principale glacier retained a surface elevation similar to that of 2010. This is visible in Fig. 5a, where it appears as an oblique line crossing the glacier from the southwest to northeast. This feature is not natural, being related to the artificial redistribution of snow to maintain ski routes on the glacier.

The average annual geodetic mass balance rate for the Dolomite glaciers varied from $-\mathrm{0.45}\pm \mathrm{0.08}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (water equivalent) in the 1980s to $-\mathrm{0.65}\pm \mathrm{0.08}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ between the 1990s and 2010 and $-\mathrm{0.63}\pm \mathrm{0.05}$ $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ in the last 13 years (Fig. 6 and Table S6).

Figure 6 Mass balance (MB) rates of World Glacier Monitoring Service (WGMS, 2021) reference glaciers (RGs) with in situ data starting in 1980 and the Dolomite glaciers. For our MB data, measurements starting in 1982 were corrected to 1980 using the average MB loss. Dolomite glacier MB with measurements starting in the 1990s were calculated back to the 1980s using the average MB rate of all Dolomite glaciers during the 1980s (i.e. 0.45 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$). Three RGs lack 2023 data (i.e. Argentière, Saint Sorlin, Serennes).

Total volume loss in the entire observation period was approximately 0.105 Gt, with 0.083 Gt lost between the 1980s and 2010 and 0.022 Gt lost during the 2010–2023 interval (Table S7). The Marmolada glacier by itself accounted for 65.7 % of total volume losses. Other major contributors included the Antelao Superiore and Fradusta glaciers with 7.9 % and 7.3 %, respectively.

Our results show that the use of total area among each comparison in the geodetic mass balance leads to an underestimation of the loss ($\mathrm{m}\mathrm{w}.\mathrm{e}.$) when compared to calculations based on common area. In our case the bias introduced by total area led to an underestimation between 1 % and 31 % of the common area geodetic mass balance, depending on the site and considered time step. There were some cases of decadal comparison (1980s–2010 in Cristallo, Antelao Inferiore, and Marmolada) where the total glacier area produced larger geodetic mass balance losses than calculations using the common area.

4.3 Climate data

The general trend for temperature and precipitation SAI showed more frequent positive values towards the end of the 1985–2022 time frame (Fig. 7a), meaning long-term above-average temperatures and precipitation. On the contrary, a negative SAI means long-term below-average temperature or precipitation. The nine lowest T SAI values fell within the first 15 years (1985–2000), with the only exception of 2014 (T SAI $=-\mathrm{1.23}$). Among the 10 highest events, 7 were in the last 15 years (2007–2022). Exceptions were the years 1994, 1998, and 2003, with the latter showing the highest T SAI with a value of 2.74. The minimum T SAI was recorded in 1989 with −2.18. Concerning the Pr SAI, the rising trend was less pronounced than for temperature but still showed a clear increase. Indeed p values were $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{2}}$ for Pr SAI and $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{5}}$ for T SAI. In this case, 6 of the 10 lowest Pr SAI values were before 2000 and 8 of the highest were recorded after 2008. The Pr SAI was greatest in 2014 with a value of >2 and lowest in 1996 with −1.22. In the first half of the time series the T SAI tended to be lower than the Pr SAI, while in the second half, this trend progressively reversed. Furthermore, SAI values, particularly regarding temperature, were more positive towards the end of the time frame. Temperatures rose by 0.4–0.6 °C per decade since 1985, while precipitation showed an increase that lasted about 15 years from 1995, culminating in the extremely snowy year of 2014 (Fig. 7b). Fedaia station, the only one providing data since 1980, did not show any trend for the total snow accumulation (p value of 0.61; Fig. 7c); however, increased extreme events could be observed in the last decade of its time frame. The snow monitoring stations did not show any trend for the total snow accumulation (p values of 0.54–0.95; Fig. 7c); however, extreme events (above the 95th percentile) were observed in 2013 and 2014 for all the stations.

Figure 7 (a) T SAI (red boxes) and Pr SAI (blue boxed) calculated between 1985 and 2022 and averaged over 22 AWSs. The boxes show the quartile (Q) of the dataset, with the central black line representing the median and green line the mean. The whiskers extend to $Q\mathrm{1}-\mathrm{1.5}\cdot (Q\mathrm{3}-Q\mathrm{1})$ and $Q\mathrm{3}+\mathrm{1.5}\cdot (Q\mathrm{3}-Q\mathrm{1})$, where Q1 and Q3 are the first and third quartiles. (b) Snow depth on the ground between October and June for the most relevant years of the study and the 1980s/1990s–2023 mean ±1σ (red line and grey shade) for the Fedaia, Ra Vales, Monti Alti di Ornella, and Col dei Baldi stations. (c) Total annual snow accumulation and 5-year centred running mean (black line) for the same snow monitoring stations. The years shown in the plot refer to hydrological years; e.g. 2023 refers to 2022/23.

Snow height on the ground showed a peak in March or April depending on the station altitude (Fig. 7b). Among the investigated years, 2023 and 1992 values were below the mean for most of the accumulation season, getting close to or even overtaking the mean only in April (April and May for the highest stations of Ra Vales and Monti Alti di Ornella).

5 Discussion

5.1 Climate trends impacting the Dolomite cryosphere

There are numerous studies describing the characteristics of the recent climate regime change occurring in the Alps (Huss et al., 2017; Hock et al., 2019). The main parameters influencing the response of the Alpine cryosphere are linked to the increase in summer temperatures and modifications in precipitation and snowfall patterns (Žebre et al., 2021). In the Eastern Alps mean summer temperature has increased by about 2.0 °C since 1979, leading to an increased ablation with more effective results also due to the extension of the ablation season (Colucci and Guglielmin, 2015; Colucci et al., 2021). Mean annual air temperature increased by 0.3±0.2 °C per decade. On the other hand, while the fraction of precipitation falling as rain has increased compared to that of snow in low- and mid-altitude areas, this effect has not yet been observed in winter in high-altitude glacial basins. Additionally, there is generally an average increase in snowfall amounts during the winter season at high altitudes in the Alps from 1971 (Matiu et al., 2021). However, the increase in snow, particularly linked to extreme events (e.g. the winter of 2014), has not been able to counterbalance the increase in summer temperatures, resulting in consistently negative mass balances after 2001 (Hugonnet et al., 2021). Only in very limited areas where the highest mean annual precipitation is recorded have some small residual glacial bodies slowed their reduction (Scotti et al., 2014; Colucci et al., 2021).

Using the SAI, it is possible to observe how in the last 2 decades unfavourable conditions for glaciation have prevailed. Only a few exceptions exist, notably 2014 which shows a winter precipitation anomaly with Pr SAI >2. That year was very favourable for glaciation in all of the Eastern Alps (Colucci et al., 2021), although the overall Alpine mass balance remained negative (WGMS, 2021). On the contrary, the mass loss of the glaciers was particularly noteworthy in the hydrological year 2021/22 (Voordendag et al., 2023). 2022 was by far the most extreme melting year among those analysed, with an extremely effective ablation season (T SAI =2.24) and a dry accumulation season (Pr SAI $=-\mathrm{1.58}$). Within Alpine mass balance records, the results of the ablation season of 2022 were unprecedented.

Overall, in the Dolomites we can observe an extreme form of the ablation period with higher mean and extreme temperatures associated with a slight increase in snowfall on glacial basins, which is however insufficient to counterbalance the temperature increase. According to such a climatic evolution, the Dolomites are rapidly turning from being mountains hosting sites favourable to local glaciation to areas where periglacial processes will progressively gain importance.

5.2 The Dolomite glaciers within the European Alps

The European Alps are home to 12 reference glaciers (RGs) with available glaciological mass balance data from before 1980 (WGMS, 2021). These can be used as benchmarks for our results, considering that the study area lacks any long-term mass balance monitoring programme. All the Dolomite glaciers exhibit a geodetic mass balance within the range of the RGs, with values at the upper end, similar to the RGs that have had less loss ($\mathrm{m}\mathrm{w}.\mathrm{e}.$) during the last 4 decades (Fig. 6). The Fradusta glacier is the only site that presents a cumulative mass balance close to the 1980–2023 RG average. Glaciers with debris cover (Popera Alto and Pensile and Sorapiss Occidentale) or which are topographically shadowed (Cristallo, Antelao Inferiore) display mass balances that are comparable to those of other RGs (e.g. Kesselwandferner, Allalin, Giétro), although these are larger and lie at higher altitudes. In terms of decadal mass balance rates, the Dolomites match the behaviour of RGs during the 1980s (−0.45 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$), while they are always lower than the average rates during the following decades. Dolomite glacier mass balance rates are half of the average RG rate during the last 13 years, with −0.63 compared to −1.27 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$. If we consider all 17 RGs available in central Europe from 2010 to 2023 the average mass balance rate equals −1.44 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (WGMS, 2021).

Calculations of glacier mass change from space are constantly evolving and improving in terms of spatial and temporal resolution (Berthier et al., 2023) and are now able to provide a valid solution for monitoring glaciers over time (Berthier et al., 2024). However, global studies with resolutions generally on the order of tens of metres (e.g. Hugonnet et al., 2021) are not suitable for monitoring very small glaciers. These glaciers, such as those in the Dolomites, have a lower degree of accuracy mainly due to their size and often present artefacts related to the surrounding steep topography. It is therefore important that global studies, which are of fundamental importance, do not replace higher-resolution measurements (e.g. close-range airborne remote sensing) for monitoring smaller glaciers.

The apparent resistance of the Dolomite glaciers to climate change is due to several factors: (i) the orographic protection offered by the complex topography of these mountains, (ii) the importance of avalanches in the dynamics of these glaciers, and (iii) the abundant debris cover that affects the remaining glacial volumes. All these factors contribute to dampening the effects of climate change, as observed in similar contexts in the Alps. Glaciers behaving similarly are for example found in the easternmost sector of the Alps (Julian Alps) or in the Orobic Alps (southern-central Alps), where topography feedback and conditions similar to the ones characterising the Dolomite glaciers are found. The Montasio Occidentale glacier and Canin East ice patches, in the Julian Alps, share less negative mass balance rates, with −0.09 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (2006–2019; De Marco et al., 2020) and −0.10 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (2010–2023; Colucci et al., 2021). The Lupo glacier, the only glacier in the Orobic Alps which is routinely monitored by the Servizio Glaciologico Lombardo (SGL), instead exhibits a mass balance rate of −1.08 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ (2010–2023), although it is heavily negatively influenced by the 2022 and 2023 seasons (Hagg et al., 2017).

Cumulative geodetic mass balance measurements derived from single years during different decades, as used in this study, can be impacted by seasonal variability influencing the surveys. The presence of residual snow on glacier surface and firn can alter the measurements introducing a bias. This bias is already partially included in the ±60 kg m⁻³ used as density conversion factor uncertainty (δ_ρ). From our observations, this might have influenced data in 1980s and the ensuing geodetic results with a negative bias in the mass balance, albeit within the methodology error. 1991–1992 glacier conditions at the time of the survey instead appear optimal. The same issues that affected the 1980s surveys may have influenced the 2010 surveys too, although with less intensity. Data from 2014 cannot be considered useful for assessing the state of glaciers, as thick residual snow was present on almost all the considered sites. Surveys from 2023 were carried out during the best possible time window and present no bias.

5.3 Local glaciological variability within the Dolomites

The glaciers of the Dolomites are all classified as very small glaciers (<0.5 km²), except for Marmolada Principale. Their short-term mass balance is influenced by local topography and avalanche accumulation, rather than by climate conditions alone.

Figure 8 Sorapiss Occidentale (a) surface elevation change and (b) 3D model. The black arrow is used as a reference between subplots (a) and (b). (c) Surface elevation change in the Marmolada glacier, with (d) a magnification of the area of interest regarding the ice avalanche that happened on 3 July 2022 with a vertical profile in 2010 and 2023 and (e) an aerial photo of the release zone.

Figure 9 Photos of different settings of the Dolomite glaciers and their relationship with local topography: (a) Popera Pensile, (b) Fradusta Superiore (multiple photos merged together), (c) frontal part of Sorapiss Occidentale, and (d) Travignolo glacier front.

An example of the influence of local topography on glacier evolution is Sorapiss Occidentale, where a steep rock wall separates the upper cirque from the lower portion of the glacier. The upper portion, where the highest thickness losses are observed, was completely debris free until a recent rockfall, which occurred after 2010 (Fig. 8a). The lower part of the glacier, where a thick layer of debris is present, extends almost up to the frontal moraine. The overload caused by the recent rockfall appears to have produced a small surge in its central crevassed portion, as proven by the advanced position of the crevasses between the 2010 and 2023 surveys. The increases in the surface elevation are related to the volume of the rock and debris fallen but cannot explain the relative elevation gain values larger than 10 m by itself (Fig. 8a). Similar processes, involving actual surges and not only partial advances, have been observed for other Alpine glaciers before on rare occasions (e.g. Deline, 2001).

The Marmolada glacier, in its main sector, Marmolada Principale, is by far the largest of the area. Thickness losses here appear to be inversely proportional to their elevation during all 4 decades of our study, although the maximum losses measured in the recent period (2010–2023) are linked to the ice avalanche detached from one of its highest portions (Fig. 8b). The former Marmolada glacier is now divided into six different parts, each with varying degrees of thickness loss (Fig. 5b). The presence of cold ice has been detected in some of them (Forte et al., 2020). It should be noted that the thickness losses of Marmolada Principale may have been reduced by anthropic activities related to the ski slope maintenance and the use of geotextile sheets, as reported during some ablation seasons (e.g. Baroni et al., 2015).

Typically, relative uncertainties of the glacier area are increasingly higher for smaller glaciers (Fischer et al., 2014). Lower uncertainties can be obtained with higher-resolution imagery (e.g. ±2.3 %, Linsbauer et al., 2021), even for regional glacier inventories. In our study, area uncertainties are not considered a relevant factor for the bare ice surface from 2010 to 2023, as we have very high-resolution multi-source datasets (i.e. point clouds, multi-temporal change, orthophotos) that do not have cloud or shadow problems and are accompanied by field observations. Therefore, the influence of area uncertainty in the subsequent mass balance calculation is limited. The issue is instead relevant to the older imagery (1980s–1990s) and for debris-covered glacier portions like Popera Alto and Sorapiss Occidentale, where mapping is not straightforward. In particular, dead-ice patches in front of the actual debris-covered glacier are not always easy to distinguish from the main glacial body if additional data (e.g. ground-penetrating radar) are not available (Santin et al., 2023) (e.g. Fig. 9a and c). Another challenging issue regards evaluating the dynamic behaviour of these debris-covered fronts. To address this, we have used the approach of Paul et al. (2009), including adjacent debris-covered areas and dead-ice bodies in the so-called “potential” glacier area. An additional verification has been done using multi-temporal data and field observations.

Some of the Dolomite glaciers, such as Fradusta Superiore, Popera Pensile, and Antelao Inferiore, are now at the border between the classification of mountain glaciers and ice patches. Despite this, they can still be classified as glaciers due to the evidence of recent dynamics and internal deformation through the presence of open crevasses. For example, Antelao Inferiore still presents visible crevasses and Popera Pensile and Fradusta Superiore show the presence of a bergschrund and annual layering (Fig. 9a and b). The outcropping ice/firn surface is often characterised by irregularly spaced convex bands of ice rich in fine sediment. These features are consistent with variations in ice/firn velocity across the glacier, with the fastest flow being observed away from the ice margins (Hughes, 2007; Colucci and Guglielmin, 2015), suggesting the presence of ice movement.

Projecting our mass balance trends into the future is not straightforward due to the high range of morphological differences in each ice body considered here. Processes like debris cover increase or glacier segmentation may contribute to different responses to ongoing warming, making these glaciers more resilient. The shrinkage and rapid disappearance of glaciers in the Dolomites may find a counterpart behaviour exhibited since the Little Ice Age in glaciers of the Julian Alps (e.g. Canin and Triglav glaciers, Colucci, 2016; Colucci and Žebre, 2016) or in other minor glacial bodies in the Prokletije (Albania), Durmitor (Montenegro), and Pirin (Bulgaria) mountain ranges (Hughes, 2010, 2018; Hughes and Woodward, 2016; Gachev et al., 2016). A similar evolution awaits all glaciers sharing a low elevation and small size in the European Alps (Cook et al., 2023).

6 Conclusions

We quantified the last 40 years of ice melting in the remaining Dolomite glaciers, highlighting how local topography can act as positive feedback in slowing mass loss compared to the larger reference glaciers of the European Alps. The existing glaciers in the region lost 56 % of their area and 0.64±0.06 $\mathrm{m}\mathrm{w}.\mathrm{e}.{\mathrm{yr}}^{-\mathrm{1}}$ from the early 1980s to 2023. The cumulative estimated volume loss is 0.105 Gt, of which 0.022 Gt was in the last 13 years (2010–2023). In the late 1950s the Dolomites hosted 33 glaciers, of which only 9 are still considered mountain glaciers today (Smiraglia et al., 2015).

All the glaciers in the Dolomites currently have their accumulation areas below the envELA and are projected to disappear or segment into minor glacial bodies in a few decades. Their fate appears inevitable even if a 1991–2020 steady-state climate is assumed. Only locally may favourable feedback be brought by the debris coverage. A few glacial bodies may eventually shift from glacial to periglacial, thus becoming more resilient in a warming climate.

Although this work adds new quantitative information that fills part of the existing gap, we strongly encourage more specific monitoring for some if not all these sites at least on an annual basis. This is of even greater importance in the context of the present climate trend, which is leading to a rapid degradation of the Dolomite glacial bodies. As evidenced by e.g. the ice avalanche of the Marmolada glacier in July 2022, the final stages of the disappearance of these glaciers may result in parossistic events, which in turn pose risks in an extremely popular tourist area.