Analysis of Regional Changes in Geodetic Mass Balance for All CaucasusGlaciers over the Past Two Decades

: Glaciers and snow in the Caucasus are major sources of runoff for populated places in many parts of this mountain region. These glaciers have shown a continuous area decrease; however, the magnitude of mass balance changes at the regional scale need to be further investigated. Here, we analyzed regional changes in surface elevation (or thickness) and geodetic mass balance for 1861 glaciers (1186.1 ± 53.3 km 2 ) between 2000 and 2019 from recently published dataset and outlines of the Caucasus glacier inventory. We used a debris-covered glacier dataset to compare the changes between debris-free and debris-covered glaciers. We also used 30 m resolution ASTER GDEM (2011) to determine topographic details, such as aspect, slope, and elevation distribution of glaciers. Results indicate that the mean rate of glacier mass loss has accelerated from 0.42 ± 0.61 m of water equivalent per year (m w.e. a − 1 ) over 2000–2010, to 0.64 ± 0.66 m w.e. a − 1 over 2010–2019. This was 0.53 ± 0.38 m w.e. a − 1 in 2000–2019. Mass loss rates differ between the western, central, and eastern Greater Caucasus, indicating the highest mean annual mass loss in the western section (0.65 ± 0.43 m w.e. a − 1 ) in 2000–2019 and much lower in the central (0.48 ± 0.35 m w.e. a − 1 ) and eastern (0.38 ± 0.37 m w.e. a − 1 ) sections. No difference was found between the northern and southern slopes over the last twenty years corresponding 0.53 ± 0.38 m w.e. a − 1 . The observed decrease in mean annual geodetic mass balance is higher on debris-covered glaciers (0.66 ± 0.17 m w.e. a − 1 ) than those on debris-free glaciers (0.49 ± 0.15 m w.e. a − 1 ) between 2000 and 2019. Thickness change values in 2010–2019 were 1.5 times more negative (0.75 ± 0.70 m a − 1 ) than those in 2000–2010 (0.50 ± 0.67 m a − 1 ) in the entire region, suggesting an acceleration of ice thinning starting in 2010. A signiﬁcant positive trend of May-September air temperatures at two selected meteorological stations (Terskol and Mestia) along with a negative trend of October-April precipitation might be responsible for the negative mass balances and thinning for all Caucasus glaciers over the study period. These results provide insight into the change processes of regional glaciers, which is key information to improve glaciological and hydrological projections in the Caucasus region.


Introduction
Glacier mass balance is an important variable for understanding the response of glaciers to climate change [1] and their contribution to regional water resources [2] as well as from the perspective of changes in global sea level [3][4][5].Glacier mass balance can be derived by two different methods, such as the traditional glaciological method, and the geodetic method.Glaciological mass balance is based on in-situ measurements of accumulation and ablation using snow pits and stakes [1], while the geodetic mass balance can be obtained from observations of elevation change derived from digital elevation models (DEMs) computed from stereo satellite images.Ground-based measurement gives glaciological mass balance relatively high accuracy at the point of measurement, which is then interpolated to the glacier area, while the main advantage of geodetic mass balance is the wide regional coverage, although, the uncertainties of geodetic mass balance measurements depend on the satellite source and is often relatively high [6].However, close range sensing such as airborne or terrestrial laser scanning significantly reduces the uncertainty to an order of decimeters to meter [7,8].Significant improvements have been made in the last decades to improve the accuracy of the geodetic approach [9,10].Furthermore, this method has become widely used in recent years due to the possibility to monitor large glacierised areas [3,4].
Meltwater of glaciers and snow provides a water source for electricity generation, ecosystems, and communities in the Greater Caucasus.These glaciers also have a positive impact on the economy by being a major tourist attraction in this region.Several national parks on both sides of the Greater Caucasus (Georgia/Russia) welcome thousands of visitors each year.On the other hand, glacier hazards are relatively common in this region, leading to major casualties [11][12][13].Thus, estimation of glacier mass change is of great significance in this region.Investigation of glacier mass balance will also provide important information for future research on various socio-economic impacts such as water resource management, risk assessments, and tourism in this mountain region.
Attempts to organize in-situ regular observations of glacier mass balance have been made since the early 1960s on both southern (Georgia) and northern (Russia) slopes of the Greater Caucasus [14,15].The duration of those measurements was short for various reasons.Long-term glaciological observation campaigns (stake method) on the Russian side of the Greater Caucasus began in 1967 on Djankuat Glacier (43 • 11 57 N; 42 • 45 24 E) [16], which started in 1982 on Garabashi Glacier (43 • 18 19" N; 42 • 28 13 E) [17].A similar monitoring program was started in 1968 on Buba (Tbilisa) Glacier (42 • 44 57 N; 43 • 44 45 E) on the Georgian side of the Greater Caucasus [18].After the dissolution of the Soviet Union in 1990s, the monitoring program was interrupted in Georgia and no glaciological field measurements were carried out until 2007, when the mass balance monitoring programs were temporarily (2007-2010) initiated on Zopkhito Glacier (42 • 52 54 N; 43 • 25 29 E).One-year glaciological mass balance measurement was also initiated on Chalaati Glacier (43 • 7 49 N; 42 • 42 30 E) in 2011 [19].The initiation of other continuous monitoring programs was difficult mainly due to hard political and social conditions resulting in more than 30 years of gaps in glaciological mass balance data.Nowadays, continuous series of long-term observations by traditional glaciological method continue only on Djankuat and Garabashi glaciers from the Russian side of the Greater Caucasus [20].
In the last few decades, research in the Greater Caucasus focused mainly on documenting changes in glacier area [21][22][23][24], indicating a significant shrinkage in regional ice coverage.Long-term ground-based measurements from Djankuat and Garabashi glaciers also suggest predominantly negative mass balances over the last few decades from this region [17,25].Recently, geodetic mass balances have also been monitored within the local [26,27], and global [3,4] studies.The main limitations of local studies are incomplete area coverage, while the global studies always describe the Caucasus, with the Middle East making the regional analysis difficult.The vast majority of glaciers in the Middle East are very small, with mean area of <0.1 km 2 [28,29] and the climate is hot and arid (dry) [30], causing a much higher rate of glacier decline [28,29].Thus, combining these two regions causes a bias with more negative values over the Caucasus region.
In this study, we use thickness and surface mass balance data from recently published global study by Hugonnet et al. (2021) [3].This study offers a unique opportunity to analyze the behavior of individual glaciers in the Greater Caucasus as the data of mass balance and thickness changes are available for each specific glacier.Thus, the main goal of this study is to: (i) analyze the behavior of glaciers in the Greater Caucasus separately from the Middle East; (ii) estimate changes in surface elevation and geodetic mass balance for all Caucasus glaciers over the past two decades according to river basins and analyze these changes on sub-regional (e.g., western vs. eastern; northern vs. southern) and individual glaciers (e.g., debris-covered vs. debris-free); (iii) compare this data with those from local existing ground-based and geodetic mass balance measurements along with local climate observations from high-mountain meteorological stations.

Study Area
Caucasus is a complex mountain system located between the Black and Caspian seas on the border of Europe and Asia consisting of two separate mountain systems such as the Greater Caucasus and the Lesser Caucasus (Transcaucasian Highlands).The Greater Caucasus is the higher and more extensive part of this region, extending for about 1100 km from northwest to southeast.The Lesser Caucasus is located approximately 100 km to the south and characterized by relatively lower elevations and warm climate, providing unfavourable conditions for the existence of modern glaciers.Thus, all contemporary glaciers of the Caucasian region belong to the Greater Caucasus.
According to geomorphological characteristics, altitudinal difference, and climatic conditions, the Greater Caucasus divided into three parts-the Western, Central, and Eastern.The central Greater Caucasus is situated between the mountain peaks of Elbrus (5642 m) and Kazbegi (5047 m) (Figure 1).Several summits in the central Greater Caucasus exceed 5000 m above sea level (a.s.l.): Dykh-Tau (5205 m), Shkhara (5203 m), Koshtan-Tau (5151 m), Jangha (Dzhangi-Tau) (5058 m), and Pushkin Peak (5034 m).The western section is the lowest part of the Greater Caucasus, with highest summit being Dombai-Ulgen (4046 m).The Eastern section of the Greater Caucasus is the most expanded part of this mountain region, with the highest summit being Mt.Bazardüzü (4466 m a.s.l.).The division of the Greater Caucasus into northern and southern macroslopes by location relative to the main crest of the range is also often used.About 13.4% of the surface area of Caucasus glaciers was covered by debris in 2014 [38].Nowadays, there are more than 2200 glaciers in the Greater Caucasus, covering an area of about 1060 km 2 .The northern slopes of the mountain range contain more glaciers than the southern slopes.About 70% of the modern glaciated area is concentrated in the central section.Fourteen glaciers have an area > 10 km 2 resulting in ~220 km 2 or ~21% of the total glacier area in 2020.The Bezingi Glacier, with an area of ~39 km 2 , was the largest The mountain range of the Greater Caucasus represents an orographic barrier to the flow of moisture brought from the west and southwest by Mediterranean and Atlantic cyclones.The greatest amount of precipitation falls on the southern slope of the western section.Here, the amount of annual precipitation at an altitude of 3000 m a.s.l.reaches 4000 mm.Precipitation gradually decreases from west to east, reaching about 1000 mm on the southern slopes of eastern Greater Caucasus.The northern slope is relatively more arid (dry) and cold than the southern one, resulting in ~500-1000 mm less annual precipitation and 1-2 • C lower summer temperatures [31,32].The average regional lapse rate is minimum in winter (−2.3 • C per 1000 m) and maximum in summer (−5.2 • C per 1000 m) [33].Overall, snow depth may reach 5-7 m in winter in several regions of the western part of the Greater Caucasus, such as the Upper Svaneti and northern Abkhazeti regions in Georgia [34].
Glacier decrease started from the last advances of the Little Ice Age (LIA) in the northern Greater Caucasus, dated approximatively to the late 1840s [35], though it may have been slightly earlier, in the 1810s, in the southern slopes [36].The decrease rate was relatively small in centennial scale (~0.4% yr −1 in 1911-2014) [22][23][24], while this has increased since the 1960s (0.5% yr −1 ) [21] and reached its maximum rates over the past two decades (1.16% yr −1 ) [37].
About 13.4% of the surface area of Caucasus glaciers was covered by debris in 2014 [38].Nowadays, there are more than 2200 glaciers in the Greater Caucasus, covering an area of about 1060 km 2 .The northern slopes of the mountain range contain more glaciers than the southern slopes.About 70% of the modern glaciated area is concentrated in the central section.Fourteen glaciers have an area > 10 km 2 resulting in ~220 km 2 or ~21% of the total glacier area in 2020.The Bezingi Glacier, with an area of ~39 km 2 , was the largest glacier mapped in the 2020 database [37].

Data and Methods
In this study, we used geodetic mass balance and ice thickness data from Hugonnet et al. (2021), [3] obtained from 30 m resolution Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) and ASTER L1A stereo images (each covering 60 km × 60 km) between 1 January 2000 and 31 December 2019 (downloaded from https://lpdaac.usgs.gov/products/ast_l1av003/)(Last access 13 December 2021).The specifically developed statistical methods are consistently described in Hugonnet et al. (2021) [3].Outlines from the Caucasus glacier inventory from 2014 [21] were regarded as glacier boundaries in Hugonnet et al. (2021) [3].These outlines for most glaciers in this area are identified according to the Landsat 8 Operational Land Imager (OLI) images and downloaded from the EarthExplorer website (http://earthexplorer.usgs.gov/)(Last access 13 December 2021).
The 5-year uncertainties by Hugonet et al., 2021 [3], at the 95% confidence level, were lower than 0.5 m yr −1 for glaciers larger than 1 km 2 and conservative for smaller glaciers.They, thus, validated the reliability of improved uncertainty approaches for volume change estimation down to the scale of individual glaciers.A detailed description of uncertainty calculation is given in Huggonet et al., 2021 [3].To check the accuracy of geodetic mass balance data by Hugonnet et al. (2021) [3] relative to the existing groundbased observation from the same period, we selected the Djankuat Glacier (Figure 2a) and extracted its glaciological mass balance data from the WGMS server [20].For more comparison with other geodetic mass balance studies, we also selected the Kyukyurtyu, Ulluchiran, Karachaul, and Bolshoy Azau glaciers from the Elbrus Massif (Figure 2b), previously studied by Kutuzov et al. (2019) [27], based on two different DEMs from 1997 (aerial photography survey) and 2017 (Pléiades).We selected only these glaciers because they have almost similar internal boundaries in Hugonnet et al. (2021) [3] and Kutuzov et al. (2019) [27] and are comparable in terms of data analysis.The rest of the glaciers from Elbrus Massif were excluded for similar comparison due to the much different internal boundaries, although the entire Elbrus Massif was also selected as a single outline with only outer/external boundaries.kyurtyu, Ulluchiran, Karachaul, and Bolshoy Azau glaciers from the Elbrus Massif (Figure 2b), previously studied by Kutuzov et al. (2019) [27], based on two different DEMs from 1997 (aerial photography survey) and 2017 (Pléiades).We selected only these glaciers because they have almost similar internal boundaries in Hugonnet et al. (2021) [3] and Kutuzov et al. (2019) [27] and are comparable in terms of data analysis.The rest of the glaciers from Elbrus Massif were excluded for similar comparison due to the much different internal boundaries, although the entire Elbrus Massif was also selected as a single outline with only outer/external boundaries.We used a debris-covered glacier dataset from Tielidze et al. (2020a) [38] to compare geodetic mass balance and thickness changes between debris-free and debris-covered glaciers.This dataset was compiled based on the Glacier Classification Guidance from the Global Land Ice Measurements from Space (GLIMS) initiative for remote-sensing observations [39].According to this guideline, all glaciers with a debris coverage of more than 10% were considered to be debris-covered and otherwise to be debris-free glaciers.We selected 34 glaciers from this dataset with an area of >2.0 km 2 .
We also used the ASTER Global Digital Elevation Model (GDEM, 17 November 2011) version 3 to determine topographic details such as aspect, slope, and elevation distribution of glaciers.The DEM was downloaded from NASA LP DAAC Collections (http://earthexplorer.usgs.gov/)(Last access 13/12/2021).
Change in the mass (M) for all Caucasus glaciers over the last 20 years was calculated based on geodetic mass balance data by Hugonnet et al. (2021) [3] and the recently published Caucasus glacier inventory by Tielidze et al. (2021) [37] based on two different satellite sources from 2000 (Landsat) and 2020 (Sentinel).
This was calculated as follows: where dm is the average surface balance (−0.53 m w.e. a  [37]. We converted an ice volume (km 3 ) to a mass of ice (in Gt) by using the following equation: Mass of ice (Gt) = Volume of ice (km 3 ) × Density of ice (Gt/km 3 or 850 ± 60 kg/m 3 ).We used a debris-covered glacier dataset from Tielidze et al. (2020a) [38] to compare geodetic mass balance and thickness changes between debris-free and debris-covered glaciers.This dataset was compiled based on the Glacier Classification Guidance from the Global Land Ice Measurements from Space (GLIMS) initiative for remote-sensing observations [39].According to this guideline, all glaciers with a debris coverage of more than 10% were considered to be debris-covered and otherwise to be debris-free glaciers.We selected 34 glaciers from this dataset with an area of >2.0 km 2 .
We also used the ASTER Global Digital Elevation Model (GDEM, 17 November 2011) version 3 to determine topographic details such as aspect, slope, and elevation distribution of glaciers.
The DEM was downloaded from NASA LP DAAC Collections (http://earthexplorer.usgs.gov/)(Last access 13 December 2021).This was calculated as follows: where dm is the average surface balance (−0.53 m w.e. a −1 ) (meter water equivalent per year, no mentioned in the rest part of this paper) for the entire Caucasus in 2000-2019 by Hugonnet et al. (2021) [3]; ∆t is the duration of the observation period (20 years); S1 and S2 are the glacier areas at the beginning (1381.5 km 2 in 2000) and at the end (1069.1 km 2 in 2019/2020) of the observation period, according to Tielidze et al. (2021) [37].
We converted an ice volume (km 3 ) to a mass of ice (in Gt) by using the following equation: Mass of ice (Gt) = Volume of ice (km 3 ) × Density of ice (Gt/km 3 or 850 ± 60 kg/m 3 ) Monthly climate data (temperature and precipitation) from two different meteorological stations were also used in this study.The Terskol station is located at an elevation of 2143 m a.s.l. on the Russian side of the Greater Caucasus, while the Mestia station is located at an elevation of 1441 m a.s.l. on the Georgian side (see location on Figure 1).Observations are available for the period of 2000-2018 for Terskol (temperatures and precipitation) and 2000-2015 for the Mestia station.

Mass Balance and Thickness Change for the Entire Region
Changes in annual geodetic mass balance along with cumulative mass balance in the entire Greater Caucasus over the last twenty years are shown in Figure 3 (see also Table A1).We computed mean annual mass loss of 0.53 ± 0.38 m w.e. a −1 for all glaciers in the Greater Caucasus between 2000 and 2019 with a cumulative mass change of −10.07 ± 7.22 m w.e.Calculated decrease in the mass of all glaciers in the Greater Caucasus over the past 20 years amounted to the ~12.99 km 3 or 11.04 ± 0.78 Gt water (~0.58Gt per year).The mean annual geodetic mass balance between 2000 and 2010 was −0.42 ± 0.61 m w.e. a −1 for the entire mountain range, while it was −0.64 ± 0.66 m w.e. a −1 between 2010 and 2019 (Figure 4).The rates of mass loss vary between the western, central, and eastern Greater Caucasus.This was relatively higher in the western section with 0.65 ± 0.43 m w.e. a −1 in 2000-2019 (Figures 4 and 5a) and much lower in the eastern section (0.38 ± 0.37 m w.e. a −1 ) (Figures 4 and 5c).The colour-coded map in Figure 5 shows spatial distribution of changes of surface mass for both slopes of western, central, and eastern Greater Caucasus glaciers in 2000-2019.No difference was found between the northern and southern slopes over the study period (0.53 ± 0.38 m w.e. a −1 on both slopes).Mean annual mass balance on the glaciers of Russian Caucasus was slightly higher (−0.55 ± 0.39 m w.e. a −1 ) than those on the Georgian counterparts (−0.51 ± 0.36 m w.e. a −1 ).Glaciers in Azerbaijan experienced the lowest mass balance (−0.23 ± 0.29 m w.e. a −1 ) over the last twenty years (Figure 4, Table A2).Mean annual surface mass balance of the north-and north-westoriented glaciers was more negative than those of the south and south-western-oriented glaciers (Figure 6a). Figure 6b shows the mean annual geodetic mass balance rates for glaciers of seven size classes of the entire mountain range between 2000 and 2019 relative to the mean and minimum (terminus) elevations, indicating that glaciers with a size class of 0.5-1.0km 2 experienced the highest mean annual mass loss (0.60 m w.e. a −1 ).The large size class of glaciers with an area of 5.0-10.0km 2 had the lowest annual mass loss (0.41 m w.e. a −1 ) during the investigated period.Changes of geodetic mass balance and thickness for individual glaciers larger than 2.0 km 2 of the Greater Caucasus in 2000-2010, 2010-2019, and 2000-2019 are shown in Table A3.

Mass Balance and Thickness Change for the Entire Region
Changes in annual geodetic mass balance along with cumulative mass balance in entire Greater Caucasus over the last twenty years are shown in Figure 3 (see also T A1).We computed mean annual mass loss of 0.53 ± 0.38 m w.e. a −1 for all glaciers in Greater Caucasus between 2000 and 2019 with a cumulative mass change of −10.07 ± m w.e.Calculated decrease in the mass of all glaciers in the Greater Caucasus ove past 20 years amounted to the ~12.99 km 3 or 11.04 ± 0.78 Gt water (~0.58Gt per year).mean annual geodetic mass balance between 2000 and 2010 was −0.42 ± 0.61 m w.e. a the entire mountain range, while it was −0.64 ± 0.66 m w.e. a −1 between 2010 and (Figure 4).The rates of mass loss vary between the western, central, and eastern Gre Caucasus.This was relatively higher in the western section with 0.65 ± 0.43 m w.e. a 2000-2019 (Figures 4 and 5a) and much lower in the eastern section (0.38 ± 0.37 m w.e (Figures 4 and 5c).The colour-coded map in Figure 5 shows spatial distributio changes of surface mass for both slopes of western, central, and eastern Greater Cauc glaciers in 2000-2019.No difference was found between the northern and sout slopes over the study period (0.53 ± 0.38 m w.e. a −1 on both slopes).Mean annual m balance on the glaciers of Russian Caucasus was slightly higher (−0.55 ± 0.39 m w.e than those on the Georgian counterparts (−0.51 ± 0.36 m w.e. a −1 ).Glaciers in Azerb experienced the lowest mass balance (−0.23 ± 0.29 m w.e. a −1 ) over the last twenty y (Figure 4, Table A2).Mean annual surface mass balance of the north north-west-oriented glaciers was more negative than those of the south south-western-oriented glaciers (Figure 6a). Figure 6b shows the mean annual geod mass balance rates for glaciers of seven size classes of the entire mountain range betw 2000 and 2019 relative to the mean and minimum (terminus) elevations, indicating glaciers with a size class of 0.5-1.0km 2 experienced the highest mean annual mass (0.60 m w.e. a −1 ).The large size class of glaciers with an area of 5.0-10.0km 2   A3.Correlation between the thinning rates and mean elevation for all Caucasus glaciers in observed periods suggests that lower elevation glaciers experienced higher annual surface lowering, i.e., the thinning rate varied significantly at an altitudinal range (Figure A1).The rate of thickness changes varies between investigated periods.We observed 1.5 times more negative values of thickness changes in 2010-2019 (0.75 ± 0.70 m a −1 ) than in 2000-2010 (0.50 ± 0.67 m a −1 ) for the entire mountain region, suggesting an acceleration of ice thinning rates starting in 2010 s (Table A1; Figure A2).The most negative values of thickness changes were observed in Russian Caucasus (0.65 ± 0.38 m a −1 ) between 2000 and 2019, while the lowest thinning rate was observed in Azerbaijan (0.27 ± 0.31 m a −1 ).Georgian glaciers experienced thickness changes of −0.60 ± 0.36 m a −1 during the same time (Table A2; Figure A2).In terms of individual river basins, the most negative elevation and mass balance change were observed between 2000 and 2010 in Bolshaya Laba basin (north-west Greater Caucasus) with 1.57 ± 1.10 m a −1 and 1.34 ± 1.20 m w.e. a −1 loss respectively (27 glaciers with a total area of 4.4 ± 0.2 km 2 ) (Table A1, Figure 1).Correlation between the thinning rates and mean elevation for all Caucasus glaciers in observed periods suggests that lower elevation glaciers experienced higher annual surface lowering, i.e., the thinning rate varied significantly at an altitudinal range (Figure A1).The rate of thickness changes varies between investigated periods.We observed 1.5 times more negative values of thickness changes in 2010-2019 (0.75 ± 0.70 m a −1 ) than in 2000-2010 (0.50 ± 0.67 m a −1 ) for the entire mountain region, suggesting an acceleration of ice thinning rates starting in 2010′s (Table A1; Figure A2).The most negative values of thickness changes were observed in Russian Caucasus (0.65 ± 0.38 m a −1 ) between 2000 and 2019, while the lowest thinning rate was observed in Azerbaijan (0.27 ± 0.31 m a −1 ).Georgian glaciers experienced thickness changes of −0.60 ± 0.36 m a −1 during the same time (Table A2; Figure A2).In terms of individual river basins, the most negative elevation and mass balance change were observed between 2000 and 2010 in Bolshaya Laba basin (north-west Greater Caucasus) with 1.57 ± 1.10 m a −1 and 1.34 ± 1.20 m w.e. a −1 loss respectively (27 glaciers with a total area of 4.4 ± 0.2 km 2 ) (Table A1, Figure 1).

Comparison with Previous Observations
Comparison with glaciological measurements from Djankuat Glacier [20] shows that values by Hugonnet et al. (2021) [3] from the period of 2007-2015 match within the uncertainty, while the rates from the other periods (2001, 2003, 2006, and 2017) show different changes (Figure 8a).In a total of twenty years, the mean (relative) difference between both methods is 55 ± 35% with a p-value of 0.172 and a standard deviation of 0.14 (Figure 8b).

Comparison with Previous Observations
Comparison with glaciological measurements from Djankuat Glacier [20] shows that values by Hugonnet et al. (2021) [3] from the period of 2007-2015 match within the uncertainty, while the rates from the other periods (2001,2003,2006, and 2017) show different changes (Figure 8a).In a total of twenty years, the mean (relative) difference between both methods is 55 ± 35% with a p-value of 0.172 and a standard deviation of 0.14   Mean annual May-September air temperatures at the Terskol weather station show a significant positive trend in the years 2000-2018, rising by 0.7 °C from 2000-2009 (10.1 °C) to 2010-2018 (10.8 °C) (Figure 9a).The warming trend is also positive for the Mestia weather station, rising by 0.9 °C from 2000-2009 (14.0 °C) to 2010-2015 (14.9 °C).A significant negative trend is seen in the October-April monthly precipitation data for the Terskol station in the years 2000-2017, decreasing by 132 mm from 2000-2009 (592 mm) to 2010-2017 (460 mm).This was also negative for the Mestia station, decreasing by 58 mm from 2000-2009 (622 mm) to 2010-2014 (564 mm) (Figure 9b).

Discussion
Results obtained from Hugonnet et al. (2021) [3] show that negative values of geodetic mass balance dominated between 2000 and 2019, but with positive values in some river basins and individual glaciers in 2000-2010 (Table A1).There was a negative and accelerated trend in the second investigated period (2010-2019) for almost all individual glaciers, with very few exceptions (e.g., Adishi, Khalde, and Kolka glaciers) (Table A3).

Discussion
Results obtained from Hugonnet et al. (2021) [3] show that negative values of geodetic mass balance dominated between 2000 and 2019, but with positive values in some river basins and individual glaciers in 2000-2010 (Table A1).There was a negative and accelerated trend in the second investigated period (2010-2019) for almost all individual glaciers, with very few exceptions (e.g., Adishi, Khalde, and Kolka glaciers) (Table A3).
Overall, observed mass loss and thickness lowering of all Caucasus glaciers with increased values since 2010 are in good agreement with recent findings of the latest glacier inventory from this region, resulting in accelerated glacier area loss over the past two decades (1.16% yr −1 ), especially since 2014 (1.85% yr −1 ) [37].

Causes of Regional Change
The decrease of surface mass balance of glaciers on different sections of the Greater Caucasus was not homogenous.The highest mean annual mass loss between 2000 and 2019 was observed in the western Greater Caucasus, while this was relatively low in the central and eastern sections.The highest mass loss in the western section can be explained by low hypsometry corresponding the lowermost mean elevation of the glaciers.A comparison of decrease rates relative to the glacier mean elevation clearly shows this linkage (Figure 10).Furthermore, western Greater Caucasus experienced the highest increase of supra-glacial debris-covered area over the past several decades [38].This part of the region is also characterized by the highest glacier retreat after the eastern Greater Caucasus [37].It is, therefore, possible that glaciers are also rapidly thinning, favouring debris-covered area over the ablation area [40,41].
The decrease of surface mass balance of glaciers on different sections of the Greater Caucasus was not homogenous.The highest mean annual mass loss between 2000 and 2019 was observed in the western Greater Caucasus, while this was relatively low in the central and eastern sections.The highest mass loss in the western section can be explained by low hypsometry corresponding the lowermost mean elevation of the glaciers.A comparison of decrease rates relative to the glacier mean elevation clearly shows this linkage (Figure 10).Furthermore, western Greater Caucasus experienced the highest increase of supra-glacial debris-covered area over the past several decades [38].This part of the region is also characterized by the highest glacier retreat after the eastern Greater Caucasus [37].It is, therefore, possible that glaciers are also rapidly thinning, favouring debris-covered area over the ablation area [40,41].Similar decrease rates between the northern and southern slopes over the study period (0.53 ± 0.38 m w.e. a −1 on both slopes) may be due to the influence of the eastern Greater Caucasus.Glaciers on the eastern Greater Caucasus only exist on the northern slopes (i.e., glaciers do not exist on the southern slopes), where we observed the lowest mass loss.This may offset (reduce) the overall regional mass balance decrease rates of the glaciers from the northern Greater Caucasus relative to their southern counterparts.
The fact that mean annual mass balance on the glaciers of the Russian Caucasus was slightly higher (−0.55 ± 0.39 m w.e. a −1 ) than those on the Georgian counterparts (−0.51 ± 0.36 m w.e. a −1 ) can be explained by more favorable conditions of Georgian glaciers due to the prevailing southwestern transport of precipitation [32].Moreover, glaciers situated in the western-central section of the Georgian Caucasus exist in high cirques with a much steeper surface.Glaciers in the western-central section of the Russian side are, on average, less steep than on the Georgian side, mainly because most valley glacier tongues in Russia are longer and reach lower altitudes than the Georgian glaciers [38], causing more ice melt on large glacier tongues on the Russian side.

Glacier Mass Change Relative to Aspect and Size Class
Glacier mass loss was in a linear relationship between aspect and mean elevation (Figure 6a), but in a non-linear relationship between size class and mean elevation.For example, glaciers with a size class of 0.5-1.0 and 1.0-5.0km 2 experienced the highest Similar decrease rates between the northern and southern slopes over the study period (0.53 ± 0.38 m w.e. a −1 on both slopes) may be due to the influence of the eastern Greater Caucasus.Glaciers on the eastern Greater Caucasus only exist on the northern slopes (i.e., glaciers do not exist on the southern slopes), where we observed the lowest mass loss.This may offset (reduce) the overall regional mass balance decrease rates of the glaciers from the northern Greater Caucasus relative to their southern counterparts.
The fact that mean annual mass balance on the glaciers of the Russian Caucasus was slightly higher (−0.55 ± 0.39 m w.e. a −1 ) than those on the Georgian counterparts (−0.51 ± 0.36 m w.e. a −1 ) can be explained by more favorable conditions of Georgian glaciers due to the prevailing southwestern transport of precipitation [32].Moreover, glaciers situated in the western-central section of the Georgian Caucasus exist in high cirques with a much steeper surface.Glaciers in the western-central section of the Russian side are, on average, less steep than on the Georgian side, mainly because most valley glacier tongues in Russia are longer and reach lower altitudes than the Georgian glaciers [38], causing more ice melt on large glacier tongues on the Russian side.

Glacier Mass Change Relative to Aspect and Size Class
Glacier mass loss was in a linear relationship between aspect and mean elevation (Figure 6a), but in a non-linear relationship between size class and mean elevation.For example, glaciers with a size class of 0.5-1.0 and 1.0-5.0km 2 experienced the highest mean annual mass loss; however they did not experience the lowest mean elevation (Figure 6b).A more linear relationship was found between size class and terminus elevation.For example, glaciers with the same size class (0.5-1.0 and 1.0-5.0km 2 ) had lower terminus elevation than the rest of the glaciers with a size class of <0.5 km 2 .The lowest terminus elevation, but not the highest mass loss, with all glaciers > 5 km 2 can be explained due to the large accumulation areas situated in high altitudes of large glaciers (e.g., Elbrus, Bezingi, Dykhsu, Karaugom, etc.).Hence, the mass loss for large glaciers occurred at the expense of glacier tongues; meanwhile, the majority of regional mass loss occurred at the expense of relatively small (<5.0 km 2 ) glaciers.

Debris-Covered and Debris-Free Glaciers
A comparison of geodetic mass balance change for debris-covered and debris-free glaciers during the investigated period showed that mean annual mass balance decrease rates are higher for debris-covered glaciers (Figure 7, Table A3).Although the debris-covered glaciers, in general, present negative rates, the supraglacial debris clearly reduces the thinning compared with debris-free surfaces at the same elevation [42].However, it is likely that altitudinal difference between those two types of glaciers can change this assumption and play a significant role in glacier thinning.For example, an average terminus elevation for the investigated debris-covered glaciers was ~200 m lower (2466 m a.s.l.) than those for the debris-free glaciers (2661 m a.s.l.), suggesting relatively warm temperatures for debris-covered glaciers based on the average regional lapse rate in summer (−0.52 • C per 100 m) [33].Supra-glacial debris thickness can also play a key role in changes of mass balance for individual glaciers.This is broadly confirmed that supra-glacial debris cover affects surface melt with increasing ablation in cases of thin debris cover (less than a few cm) or decreasing ablation under continuous thick debris cover [43][44][45].It is also possible that ice cliffs and supra-glacial ponds play a key role by enhancing the total ablation of debris-covered glaciers in the Greater Caucasus.Ice cliffs have been identified as an important contribution to the ablation of debris-covered glaciers, even when their spatial extent is very small [46].This phenomenon was studied by terrestrial photogrammetry and unmanned air vehicle (UAV) surveys confirming that ice cliffs have a net ablation rate about 3 times higher than the average glacier tongue surface [42].We note that further studies with high-resolution imagery (e.g., UAV, LiDAR) are required to understand the influence of ice cliffs in the mass balance for individual debris-covered glaciers, or at the basin and region scale in the Greater Caucasus.

Difference between the Global and Local Studies
The difference between Hugonnet et al. (2021) [3] and Kutuzov et al. (2019) [27] from Elbrus glaciers may be due to the high-resolution DEMs used by Kutuzov (2019) [27] and the internal boundaries of the selected glaciers, which were very similar, but not the same.In addition, the study by Hugonnet et al. (2021) [3] covers 2000-2019, which is slightly different (1997-2017) than in Kutuzov et al. (2019) [27].Overall, the difference between those measurements match within the uncertainties, while the difference in geodetic mass balance from Hugonnet et al. (2021) [3] and glaciological mass balance from ground-based measurement [20] of the Djankuat Glacier was relatively high.This can be due to the much higher accuracy of ground-based measurement [47] or lower accuracy of geodetic measurement [3].However, the results from those two measurements match within the uncertainties and tightly match the first time-period (2000-2010), suggesting that they are comparable in terms of long-time measurement.

Glacier Mass Balance Relative to Local Climate Change
Glacier mass balance in the Greater Caucasus is very sensitive to summer temperatures and winter precipitation [47].Accelerated mass balance decrease from the second investigated periods (2010-2019) may be explained by a sustained increase in summer temperatures and a decrease in winter precipitation associated with regional global warming.The average summer temperature in the mountainous region of the Greater Caucasus has increased by 0.5-0.7 • C over the past 30 years [48], which was also confirmed from our observation from Terskol and Mestia meteorological stations (Tables 1 and 2; Figure 11).In parallel of the increased summer temperatures, extension of melting period was observed in the northern Greater Caucasus.For example, melting on the Garabashi Glacier at the snowline altitude (3800-4000 m a.s.l.) ended up 2-3 weeks later (in early October) over the recent years [17].Furthermore, the extension of ablation season over the last two decades was also confirmed by instrumental measurement from the southern Greater Caucasus (Zopkhito and Chalaati glaciers) (Figure 7 in [19]).The increased incoming shortwave radiation, observed since the 1980s (10 W/m 2 over 10 years), may also play a triggering role in the accelerated mass loss of the Caucasus glaciers over the past several decades.It has been proposed that this trend is associated with a weakening of the processes of formation of high and low clouds, which is due to an increase in the frequency of anticyclones in the warm season [49].
Snow-ice-covered areas of the Greater Caucasus are subject to dust outbreak events originating from desert regions able to significantly decrease snow albedo.Several dust depositions events that originated from the Sahara between the 2009-2012 and 2018 period were observed in this region [50,51].Dust events systematically shortened the snow cover duration, especially in the western Greater Caucasus.The shortening was higher for locations with higher accumulation and higher elevation [51] considering one of the significant indicators for glacier thinning in this region.
The decrease in albedo, along with increased temperatures over the ablation period causes up-glacier migration of the snowline.This has a strong impact on small glaciers that have practically lost their firn reserves over the past two decades.Consequently, the firn can no longer protect the ice, and the process of ice melting becomes even more intense due to the lower albedo (Figure 12).

Conclusions
In this study, we used glacier mass balance and thickness data obtained from geodetic method by Hugonnet et al. (2021) [3] to clarify the response of region-wide glacier mass balance to climate change in the Greater Caucasus over the last two decades (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019).We also analyzed the long-term spatio-temporal glacier surface elevation and mass change for the regional, river basin, and individual glaciers.Besides, we provided a key dataset (see Supplement) that can be useful for further studies, including, for example, hydrological modeling of runoff changes, as well as for the implementation of water management plans in this region.The main limitation of this study is related to the resolution of the data.This generally affects small glaciers more than larger complexes.Applying relatively high uncertainties was a conservative solution to this issue.This can be improved by using the high-resolution stereo imagery and continuous ground-based investigation in the future.However, the estimates of Djankuat and Elbrus glaciers from this dataset agree with glaciological and geodetic mass balance measurements from previous studies.
The main results indicated a mean regional glacier mass balance change of −0.53 ± 0.38 m w.e. a −1 in 2000-2019 with an equivalent cumulative mass change of −10.07 ± 7.22 m w.e.Calculated loss in the mass of all glaciers in the Greater Caucasus over the past 20 years amounted to the ~12.99 km 3 or 11.04 ± 0.78 Gt water.Strong acceleration in mass loss and the thinning rates were observed since 2010, which likely occurred due to an increase in summer temperatures, accompanied by a reduction in winter precipitation.
The highest mass loss in the western Greater Caucasus was due to the lowest mean elevation of glaciers.We also hypothesized that the deposition events that originated from the Sahara between the 2009-2012 and 2018 periods may have increased glacier melting in the western section.
Lower terminus elevation in debris-covered glaciers indicated a strong decrease of surface mass balance relative to the high-terminated debris-free glaciers.Debris thickness, ice cliffs, and supraglacial ponds may also play a key role in enhancing the total ablation of debris-covered glaciers in the Greater Caucasus.) Elevation (m a.s.l.)

Figure 1 .
Figure 1.Glacier extent in 2014 in the Greater Caucasus according to the main river basins.Green circles correspond to changes in glacier surface mass balance (m w.e. a −1 ) according to individual river basins.See also Table A1 to compare mass balance data relative to all river basins.

Figure 1 .
Figure 1.Glacier extent in 2014 in the Greater Caucasus according to the main river basins.Green circles correspond to changes in glacier surface mass balance (m w.e. a −1 ) according to individual river basins.See also Table A1 to compare mass balance data relative to all river basins.

Figure 2 .
Figure 2. Location of Djankuat (a) and Elbrus (b) glaciers.Glaciers with blue colour from the Elbrus Massif were used for further comparison (see Chapter 4.2).SPOT image 2019 and ALOS PALSAR DEM was used as a background.See also Figure 1 for relative location of these glaciers to the Greater Caucasus.

Figure 2 .
Figure 2. Location of Djankuat (a) and Elbrus (b) glaciers.Glaciers with blue colour from the Elbrus Massif were used for further comparison (see Chapter 4.2).SPOT image 2019 and ALOS PALSAR DEM was used as a background.See also Figure 1 for relative location of these glaciers to the Greater Caucasus.
Change in the mass (M) for all Caucasus glaciers over the last 20 years was calculated based on geodetic mass balance data by Hugonnet et al. (2021) [3] and the recently published Caucasus glacier inventory by Tielidze et al. (2021) [37] based on two different satellite sources from 2000 (Landsat) and 2020 (Sentinel).
had lowest annual mass loss (0.41 m w.e. a −1 ) during the investigated period.Changes o odetic mass balance and thickness for individual glaciers larger than 2.0 km 2 of Greater Caucasus in 2000-2010, 2010-2019, and 2000-2019 are shown in Table

Figure 4 .
Figure 4. Boxplots based on data of the individual glaciers showing changes in glacier annual geodetic mass balance (m w.e. a −1 ) for all Caucasus glaciers according to the regions (western, central, and eastern), slopes (northern, southern, entire) and countries in 2000-2010 (blue), 2010-2019 (brown), and 2000-2019 (grey).Box gives lower and upper quartiles (vertical stems) along with the mean (black horizontal line inside of the box) and median (X symbol inside of the box) value of glacier mass change rate.The circles correspond to the outliers.

Figure 4 .
Figure 4. Boxplots based on data of the individual glaciers showing changes in glacier annual geodetic mass balance (m w.e. a −1 ) for all Caucasus glaciers according to the regions (western, central, and eastern), slopes (northern, southern, entire) and countries in 2000-2010 (blue), 2010-2019 (brown), and 2000-2019 (grey).Box gives lower and upper quartiles (vertical stems) along with the mean (black horizontal line inside of the box) and median (X symbol inside of the box) value of glacier mass change rate.The circles correspond to the outliers.

Figure 6 .
Figure 6.Changes in mean annual mass balance (m w.e. a −1 ) vs. mean and minimum elevation according to the aspect (a) and seven size classes (b) for all Caucasus glaciers in 2000-2019.Vertical stems show an uncertainty estimation of the mean values.

Figure 6 .
Figure 6.Changes in mean annual mass balance (m w.e. a −1 ) vs. mean and minimum elevation according to the aspect (a) and seven size classes (b) for all Caucasus glaciers in 2000-2019.Vertical stems show an uncertainty estimation of the mean values.

4. 3 .Figure 7 .
Figure 7. Changes in mean annual mass balance (m w.e. a −1 ) for debris-free (a) and debris-covered (b) glaciers vs. terminus elevation in 2000-2019.Red dots indicate the terminus elevation of each glacier.Vertical stems show an uncertainty estimation of the mean values.

Figure 7 .
Figure 7. Changes in mean annual mass balance (m w.e. a −1 ) for debris-free (a) and debris-covered (b) glaciers vs. terminus elevation in 2000-2019.Red dots indicate the terminus elevation of each glacier.Vertical stems show an uncertainty estimation of the mean values.

Figure 8 .
Figure 8.Comparison of mean annual mass balance (m w.e. a −1 ) (a) and cumulative mass balance (m w.e.)(b) for Djankuat Glacier using glaciological (blue) [20] and geodetic (brown) [3] methods.The light gray color corresponds to the uncertainty of the geodetic method.Comparison of mean annual mass changes (m w.e. a −1 ) for Elbrus glaciers using geodetic methods (c) by Kutuzov et al. (2019) [27] (blue) and Hugonnet et al. (2021) (brown) [3].Vertical stems on panel c show an uncertainty estimation of the mean values.

Figure 8 .
Figure 8.Comparison of mean annual mass balance (m w.e. a −1 ) (a) and cumulative mass balance (m w.e.)(b) for Djankuat Glacier using glaciological (blue) [20] and geodetic (brown) [3] methods.The light gray color corresponds to the uncertainty of the geodetic method.Comparison of mean annual mass changes (m w.e. a −1 ) for Elbrus glaciers using geodetic methods (c) by Kutuzov et al. (2019) [27] (blue) and Hugonnet et al. (2021) (brown) [3].Vertical stems on panel c show an uncertainty estimation of the mean values.

Figure 9 .
Figure 9. Mean annual May-September air temperatures (a) and sum of the October-April precipitation (b) at the Terskol and Mestia meteorological stations over the last twenty years.The dotted lines are the trend showing a temperature/precipitation course with time.

Figure 9 .
Figure 9. Mean annual May-September air temperatures (a) and sum of the October-April precipitation (b) at the Terskol and Mestia meteorological stations over the last twenty years.The dotted lines are the trend showing a temperature/precipitation course with time.

Figure 10 .
Figure 10.Change in mean annual geodetic mass balance (m w.e. a −1 ) as a function of mean elevation (m a.s.l.) in western (a), central (b), and eastern (c) Greater Caucasus between 2000 and 2019.Elevation bands in 200 m.Vertical stems show an uncertainty estimation of the mean values.

Figure 10 .
Figure 10.Change in mean annual geodetic mass balance (m w.e. a −1 ) as a function of mean elevation (m a.s.l.) in western (a), central (b), and eastern (c) Greater Caucasus between 2000 and 2019.Elevation bands in 200 m.Vertical stems show an uncertainty estimation of the mean values.

Figure 11 .
Figure 11.Average May-September monthly air temperatures for the Terskol meteorological station in 2000-2009 and 2010-2018 (a) and for the Mestia meteorological station in 2000-2009 and 2010-2015 (b); Sum of the October-April precipitation for the Terskol meteorological station in 2000-2009 and 2010-2017 (c) and for the Mestia meteorological station in 2000-2009 and 2010-2014 (d).

Figure 11 .
Figure 11.Average May-September monthly air temperatures for the Terskol meteorological station in 2000-2009 and 2010-2018 (a) and for the Mestia meteorological station in 2000-2009 and 2010-2015 (b); Sum of the October-April precipitation for the Terskol meteorological station in 2000-2009 and 2010-2017 (c) and for the Mestia meteorological station in 2000-2009 and 2010-2014 (d).

Figure 11 .
Figure 11.Average May-September monthly air temperatures for the Terskol meteorological station in 2000-2009 and 2010-2018 (a) and for the Mestia meteorological station in 2000-2009 and 2010-2015 (b); Sum of the October-April precipitation for the Terskol meteorological station in 2000-2009 and 2010-2017 (c) and for the Mestia meteorological station in 2000-2009 and 2010-2014 (d).

Figure A1 .Figure A2 .Figure A2 .
Figure A1.Change in in glacier thickness or surface elevation (m a −1 ) as a function of glacier mean elevation (m a.s.l.) in the Greater Caucasus between 2000 and 2019.Elevation bands in 200 m.Vertical stems show an uncertainty estimation of the mean values.Atmosphere 2022, 13, 256 22 of 24

Table 1 .
Mean monthly air temperatures ( • C) for the Terskol and Mestia meteorological stations.

Table 2 .
Sum of the monthly precipitation (mm) for the Terskol meteorological station in 2000-2009 and 2010-2017, and for the Mestia meteorological station in 2000-2009 and 2010-2014.