Strong acceleration of glacier area loss in the Greater Caucasus between 2000 and 2020

An updated glacier inventory is important for understanding glacier behaviour given the accelerating glacier retreat observed around the world. Here, we present data from a new glacier inventory for two points in time (2000, 2020) covering the entire Greater Caucasus (Georgia, Russia, and Azerbaijan). Satellite imagery (Landsat, Sentinel, SPOT) was used to conduct a remote-sensing survey of glacier change. The 30 m resolution Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM; 17 November 2011) was used to determine aspect, slope, and elevations, for all glaciers. Glacier margins were mapped manually and reveal that in 2000 the mountain range contained 2186 glaciers with a total glacier area of 1381.5 ± 58.2 km². By 2020, the area had decreased to 1060.9 ± 33.6 km² a reduction of 23.2 ± 3.8 % (320.6 ± 45.9 km²) or −1.16 % yr⁻¹ over the last 20 years in the Greater Caucasus. Of the 2223 glaciers, 14 have an area > 10 km², resulting in the 221.9 km² or 20.9 % of total glacier area in 2020. The Bezengi Glacier with an area of 39.4 ± 0.9 km² was the largest glacier mapped in the 2020 database. Glaciers between 1.0 and 5.0 km² accounted for 478.1 km² or 34.6 % in total area in 2000, while they accounted for 354.0 km² or 33.4 % in total area in 2020. The rates of area shrinkage and mean elevation vary between the northern and southern and between the western, central, and eastern Greater Caucasus. Area shrinkage is significantly stronger in the eastern Greater Caucasus (−1.82 % yr⁻¹), where most glaciers are very small. The observed increased summer temperatures and decreased winter precipitation along with increased Saharan dust deposition might be responsible for the predominantly negative mass balances of Djankuat and Garabashi glaciers with long-term measurements. Both glacier inventories are available from the Global Land Ice Measurements from Space (GLIMS) database and can be used for future studies.

1 Introduction

Glaciers are retreating and losing mass in most regions of the world, largely in response to the ongoing atmospheric warming (Hock et al., 2019; Zemp et al., 2019; Hugonnet et al., 2021). This knowledge can only be obtained when a baseline dataset (a glacier inventory) is available to calculate glacier-specific information. Complete and accurate glacier inventories also provide the information required for various hydrological and climate modelling applications (Vaughan et al., 2013) as well as change assessment. Accordingly, a frequent update of glacier inventories is required to reduce uncertainties in subsequent calculations (Paul et al., 2020). Updated glacier inventories are also critical to outline environmental policies for glacier protection and monitoring programmes, as well as for developing mitigation and adaptation strategies in response to the impact of climate changes on future glacier development (Pfeffer et al., 2014; Huss et al., 2017).

Glaciers are an important source of fresh water in countries of the Caucasus region, and runoff in large glacier-fed rivers supplies several hydroelectric power stations. They are also important reservoirs of water for the population living downstream, often providing meltwater during seasonal droughts. Furthermore, glaciers play a significant role in the economy of the Caucasus countries as a major tourist attraction with thousands of visitors each year. Finally, they are the source of or contribute to severe natural hazards in this region (complete detachment of ice and rock, glacier surging, glacier lake outburst floods) (Evans et al., 2009; Chernomorets et al., 2018; Tielidze et al., 2019), requiring a good understanding of related processes to reduce the impact of future events on human well-being. Thus, the comprehensive study of the Caucasus glaciers is crucial for the scientific study of climate change impacts but also for societal applications or sustainable regional development.

Glaciers of the Greater Caucasus started decreasing from their Little Ice Age (LIA) maximum extent in the first half of 19th century (Solomina, 2016; Tielidze et al., 2020a), reaching the highest decrease rates (∼ 0.5 % yr⁻¹) over the past decades (Shahgedanova et al., 2014; Tielidze and Wheate, 2018). A continued decrease in Caucasus glaciers could also lead to considerable changes in glacier runoff, with implications for regional water resources. Therefore, continued glacier inventorying across this region is essential. This will also potentially reduce the uncertainties for further climatic and hydrological modelling in this region as consistent multi-temporal glacier outlines are a key input for calibration and/or validation of glacier evolution models.

In this study we present two new glacier inventories (from 2000 and 2020) for the Greater Caucasus region derived from multi-temporal optical satellite images (Landsat, Sentinel-2, SPOT 6/7) in combination with digital elevation models (DEMs) along with the observed changes. We also compare the new inventories with those already available from public databases such as the Global Land Ice Measurements from Space (GLIMS) and version 6 of the Randolph Glacier Inventory (RGIv6) and highlight the related improvements.

The year 2000 inventory was compiled following the demand for creating improved glacier outlines as close as possible to that year for version 7 of the RGI. It was created because satellite images with the required quality were available from Landsat 5 and 7. The year 2020 inventory was created to also test the improved quality of the 10 m resolution Sentinel-2 data and compare results against even higher-resolution data from SPOT6/7 and Google Earth.

Figure 1 (a) The extent of modern glaciers in the Greater Caucasus. (b) Landsat 5 TM and Landsat 7 ETM+ satellite scenes from 1999–2002. (c) Sentinel-2 satellite scenes from 2019–2020. (d) SPOT 6/7 satellite scenes from 2019. The insert map in the upper right shows the location of the Caucasus region (© Google Earth).

2 Study area

2.1 General characteristics

The Greater Caucasus mountain range is situated between the Black and Caspian seas and stretches for about 1300 km from west-northwest to east-southeast. Its width ranges from 30 to 180 km. The average elevation for its western, central, and eastern sectors is 3200, 4100, and 3700 m, respectively, with the highest point being Mt. Elbrus (5642 m). The highest central sector is situated between Mt. Elbrus and Mt. Kazbegi (5047 m), with at least five other peaks exceeding 5000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1). Almost 70 % of the Caucasus glaciers are situated in the central section. About 13.4 % of the surface area of 659 glaciers was covered by debris in 2014 (Tielidze et al., 2020b).

The Greater Caucasus is in the path of the Mediterranean and Atlantic cyclones, which carry moisture from the west and southwest. The maximum amount of precipitation falls in the southern slope of the western region, with annual precipitation of about 3200 mm. This amount declines to 2000 mm in the central section and to 1000 mm in the eastern part (Volodicheva, 2002). The mean annual temperatures at the southern slopes are usually 1–2 ^∘C higher than those in the north (Tielidze and Wheate, 2018). At the mean elevation of glaciers (around 3400 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), they are around −5.0 ^∘C (Kutuzov et al., 2016; Tielidze, 2016). The average regional lapse rate has a maximum in summer (−5.2 ^∘C per km) and a minimum in winter (−2.3 ^∘C per km) (Kozachek et al., 2017).

2.2 Previous studies

The Caucasus is one of the most studied glacierized regions in the world. The first information about glaciers dates back to the 18th and 19th centuries (Kotlyakov et al., 2015; Tielidze, 2016). The first inventory of the Caucasus glaciers was published at the beginning of the 20th century (Podozerskiy, 1911). This was the result of the compilation of a topographic map, which was carried out by military topographers from 1881 to 1910 (1329 glaciers, with a total area of 1967.4 km²). Based on the same maps and in situ measurement, Reinhardt (1916) determined a summer snowline elevation at ∼ 3100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ in the Georgian Caucasus.

The second inventory of the Caucasus glaciers was initiated within the framework of the International Hydrological Decade (1965–1975) (Catalog of Glaciers of the USSR, 1967–1978; Vinogadov et al., 1978). This inventory was created based on aerial photographs from 1955–1960, topographic maps from 1960s, and data from field observations. The inventory does not contain digital outlines of glaciers but includes only tables with glacier parameters (2002 glaciers, with a total area of 1421.78 km²). Based on the same aerial imagery and topographic maps, the elevation of summer snowline increased from about 2600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (western section) to 3600 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (eastern section). On the Elbrus Massif the snowline reached 3800 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (World Atlas of Snow and Ice Resources, 1977). At the same time, the snowline for the Georgian Caucasus was measured at ∼3270 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, with the highest values (∼3470 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) in the eastern Georgian Caucasus (Gobejishvili, 1995; Tielidze, 2017).

The state of the Caucasus glaciers was determined within the framework of the Global Land Ice Measurements from Space (GLIMS) initiative using the ASTER and Landsat (1999–2004) satellite images (Khromova et al., 2016). The number and area of glaciers were calculated only for 21 river basins (out of 53), and an incomplete but first digital database was created for the southern and northern slopes of the Greater Caucasus (1706 glaciers, with a total area of 1174.52 km²). This database was later also used for version 6 of the Randolph Glacier Inventory (RGIv6) that incorporated nominal glaciers (circles covering an area equivalent to glacier size) in the eastern and western Caucasus sections (Tielidze and Wheate, 2018) from the World Glacier Inventory – Extended Format (WGI-XF; Cogley, 2009).

Recently, an updated and expanded glacier inventory covering the entire Greater Caucasus was compiled by Tielidze and Wheate (2018). The authors used large-scale topographic maps and satellite imagery (Corona, Landsat 5, Landsat 8, and ASTER) to conduct a remote-sensing survey of glacier change for three time periods (1960, 1986, 2014), with a total glacier area of 1193.2 ± 54.0 km² in 2014.

Figure 2 (a) Three glaciers shown on the Landsat scene (28 July 2000), which were later selected for multiple digitizing. An example of multiple digitizing based on the imagery from the same year – (b) Sentinel-2 23 August 2019, (c) SPOT 13 August 2019, and (d) Google Earth 14 September 2019. Image credits: (a) https://earthexplorer.usgs.gov (last access: November 2020), (b) Copernicus Sentinel data 2019, (c) Azercosmos satellite data, (d) © Google Earth 2022.

3 Data sources

We processed eight Landsat 5 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) scenes from 1999–2002 along with nine Sentinel-2 and five SPOT 6/7 scenes from 2019–2020 to cover the entire study region in both periods (Fig. 1, Table S1 in the Supplement). In addition, high-resolution QuickBird images (2019) superimposed upon the SRTM3 topography (Raup et al., 2014) were used through to the Google Earth. All the Landsat and Sentinel scenes were downloaded from EarthExplorer (http://earthexplorer.usgs.gov) (last access: November 2020), while the orthorectified high-resolution (spatial resolution 1.5 m) SPOT scenes were received from Azercosmos (https://azercosmos.az/, last access: November 2020). The Sentinel scenes served as a basis for glacier mapping, while the Google Earth and SPOT scenes were used for correction of glacier outlines and comparison of manually mapped glacier margins to those of Sentinel-2 scenes from the same year (Fig. 2) (see also Sect. 4.1). All images were acquired at the end of the ablation season, ranging from 28 July to 12 September, when glaciers were mostly free of seasonal snow under cloud-free conditions. In the case of local clouds, shadow, or snow cover, a few additional scenes from the same period were used to correctly digitize glacier outlines.

We used false-colour composites for each Landsat acquisition date, combining the shortwave infrared (SWIR), near-infrared (NIR), and red bands as RGB. The panchromatic band (15 m resolution) from Landsat 7 ETM+ was also used for better identification of glacier extents. For Sentinel-2, the colour composites were created from 10 m resolution visible and near-infrared band composites, resulting in much higher quality outlines than those derived from the Landsat scenes. The 20 m resolution SWIR band (11) was bilinearly resampled to 10 m resolution to obtain glacier outlines at this resolution automatically (e.g. Paul et al., 2020).

The ASTER Global Digital Elevation Model (GDEM, 17 November 2011) version 3 was used 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: November 2020).

Figure 3 (a) Selected glaciers for multiple digitizing based on Landsat 5 TM scene (12 August 2000). (b) Manually mapped debris-free glacier outline with a half-pixel (15 m) buffer interval. (c) Manually mapped debris-covered glacier outline with a two-pixel (60 m) buffer interval (light blue). Image credit: (a) https://earthexplorer.usgs.gov (last access: November 2020).

4 Methods

4.1 Glacier mapping

Glacier boundaries have been manually delineated from our study area. This mapping method is well adopted for the Caucasus region (e.g. Shahgedanova et al., 2014; Tielidze, 2016; Tielidze and Wheate, 2018) despite some advantages of the automated mapping method of clean ice (Paul et al., 2013). This decision was made due to the significant amount of debris-covered glaciers in this region (Tielidze et al., 2020b) as well as deep shadows where automated mapping often fails and manual corrections are required (Paul et al., 2013). Moreover, seasonal snow off glaciers was present in several scenes, and instead of removing them after an automated classification they were just not digitized. We acknowledge that identification of such non-glacier snow patches was sometimes difficult and is a highly subjective process. As a guide, we excluded snow-only features and those with a complex perimeter. This was facilitated by local experience and using the outlines of the previous glacier inventory (Tielidze and Wheate, 2018). The size of the smallest glacier included was finally restricted to 0.01 km². Glacier length was measured from changes of a centre line (Paul and Svoboda, 2009).

Estimation of the glacier mapping uncertainty is necessary to assess the significance of derived glacier changes and avoid misinterpretation of mapping. For this purpose, first we tested multiple digitization as a supplementary tool for uncertainty assessment of glacier margin identification (Paul et al., 2013). A sub-sample of three glaciers from a high-resolution SPOT image with areas of 0.3–6.3 km² were re-digitized by three different operators. The selected glaciers included Maili (42^∘43^′21^′′ N, 44^∘28^′36^′′ E), Chachi (42^∘43^′14^′′ N, 44^∘30^′20^′′ E), and G044493E42730N (GLIMS ID) (42^∘43^′41^′′ N, 44^∘29^′40^′′ E). All outlines obtained from the SPOT image (13 August 2019), along with original (basic) outlines obtained from the Sentinel image (23 August 2019), were then exported to Google Earth (14 September 2019) for comparison and visual inspection (Fig. 2b–d). The uncertainty for two debris-free glaciers (Chachi and G044493E42730N) based on normalized standard deviation (NSD – delineations by multiple digitalization divided by the mean glacier area for all outlines) was small at 1.8 % while the one debris-covered glacier (Maili) showed a much higher uncertainty at 5.1 %. The average uncertainty between the two datasets was calculated as 3.5 %. A similar approach was used for glaciers ranging from 0.4 to 6.1 km² from the Landsat imagery. The selected glaciers included Kirtisho (42^∘49^′52^′′ N, 43^∘35^′37^′′ E), Bartuy (42^∘49^′54^′′ N, 43^∘37^′33^′′ E), Khvargula (42^∘48^′12^′′ N, 43^∘37^′29^′′ E), and four relatively small neighbouring glaciers. The mapping uncertainty for debris-free glaciers was 2.1 %, while it was 6.7 % for debris-covered glaciers and 4.4 % for all glaciers of this sample (Fig. 3).

Table 1 The Greater Caucasus glacier count and area change in 2000–2020 by country.

Figure 4 Relative frequency histograms for (a) glacier area and (b) count for the seven glacier size classes in the Greater Caucasus in 2000 and 2020.

We used the buffer method as a further tool of uncertainty estimation for the entire Greater Caucasus. A buffer was drawn around the glacier outlines using ArcGIS 10.6.1 software, as suggested by Granshaw and Fountain (2006). For the images of 2020 we used a buffer equal to the resolution of the Sentinel scenes (10 m) and a half-pixel-sized buffer (15 m) for the glacier outlines derived from Landsat images for 2000. The selected buffer size for Landsat scenes is based on a recent study from the Caucasus region (Tielidze et al., 2020b) while the Sentinel buffer was selected based on a study from the European Alps (Paul et al., 2020). We assume that the larger buffer should be used for debris-covered parts of the glaciers, due to their higher uncertainty (Tielidze et al., 2020b). However, we did not enforce this here, as the related calculations are computationally difficult and challenging (Mölg et al., 2018) and would still not reflect the real problem in debris identification (see Paul et al., 2020). Instead, we used buffer with a size of two pixels for debris-covered glaciers (e.g. Frey et al., 2012), resulting in an upper-bound value of the uncertainty (Paul et al., 2020) (Fig. 3). Overall, the mapping uncertainty of the total glacier area was calculated as ±33.6 km² or (±3.2 %) for Sentinel data from 2020, which is comparable with our uncertainty estimate based on the multiple digitization method (±3.5 %). For Landsat data from 2000, the buffer uncertainty was calculated as ±58.2 km² or ±4.3 %, again comparable with the multiple digitization method for Landsat imagery (±4.4 %). It was explored that the larger glacier outlines had smaller uncertainty than the small glaciers.

Other potential uncertainties were related to the interpretation and manual digitization of the glacier margins (e.g. seasonal snow, topographic shadows, and supraglacial debris). To reduce the effect of this uncertainty, local knowledge and outlines from a previous glacier inventory (Tielidze and Wheate, 2018) were used as a delineation reference source.

4.2 Terminus measurement

Changes in the glacier terminus are a delayed and filtered response to changes in climate and are thus widely used to demonstrate climate change impacts for a large public (Lea et al., 2014). Their interpretation in climatic terms is, however, challenging as glacier-specific characteristics (e.g. response times) have to be considered (Oerlemans, 2005). Front variation measurements were conducted by intersecting the glacier outlines for each date with the centre lines. Additionally, we measured the elevations of the point at the intersection to determine the change in elevation of the glacier fronts. Length change uncertainties for the related glaciers were calculated according to source image resolution following Hall et al. (2003).

Table 2 Cumulative glacier area and count change for seven size classes in the Greater Caucasus by slope and section in 2000–2020. Bold numbers indicate the initial size class of glaciers in 2020 to fairly determine the decrease in area per size class between 2000 and 2020, while the other numbers show the absolute glacier area and count in 2020 by the same size classes (see also Fig. 9).

Figure 5 The histogram of glacier area distribution along with mean elevation (dotted line) in 2000 and 2020 for (a) the northern, (b) the southern, and (c) the entire Greater Caucasus.

Figure 6 Proportion of glacier aspect by (a) count and (b) area (km²) in 2000–2020.

5 Results

5.1 Glacier inventory 2000

Based on Landsat data from 2000 we have identified and mapped 2186 glaciers larger than 0.01 km² with a total area of 1381.5 ± 58.2 km² from 53 river basins in the Greater Caucasus (Table S2 in the Supplement). From this, 931.6 ± 37.7 km² or 67.4 % of the total glacier area was found in Russia, 446.6 ± 19.9 km² or 32.3 % in Georgia, and 3.4 ± 0.3 km² or 0.3 % in Azerbaijan (Table 1). The mean glacier size for the entire mountain region was 0.63 km², and the glacier size class 1.0–5.0 km² dominated with a total area of 478.1 km² (Fig. 4a), which is 34.6 % of the total by area. The glacier size class 0.1–0.5 km² had the most when counting by number (837 glaciers) in 2000 (Table 2, Fig. 4b). The pattern of size classes was different in the western Greater Caucasus compared to those in the central and eastern parts. The mean elevation of the glaciers ranged from 3300 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (southern slope) to 3480 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (northern slope), with an average of 3430 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 5). The number distribution by aspect showed that the glaciers were predominantly oriented towards the northwest (538 glaciers), while according to the area, the majority of the glaciers were oriented northeast (330.4 km²) (Fig. 6).

The total area of 20 glaciers from Elbrus Massif was mapped as 121.5 ± 2.2 km² in 2000 (Fig. S1 in the Supplement). The three largest glaciers mapped from the Greater Caucasus based on Landsat imagery (2000) are Bezengi – 39.4 ± 0.9 km² (43^∘2^′47^′′ N, 43^∘4^′0^′′ E), Dykhsu – 33.6 ± 0.9 km² (42^∘59^′5^′′ N, 43^∘10^′46^′′ E) (Russia), and Lekhziri 32.8 ± 0.9 km² (43^∘9^′26^′′ N, 42^∘45^′54^′′ E) (Georgia).

5.2 Glacier inventory 2020

Over the entire Greater Caucasus, the total glacier area mapped for 2020 is 1060.9 ± 33.6 km² (2223 glaciers) (Table S2). From this, 719.4 ± 22.9 km² (67.8 %) of glacier area is found in Russia, 340.8 ± 11.2 km² (32.1 %) in Georgia, and 0.8 ± 0.04 km² (0.1 %) in Azerbaijan (Table 1). Very small glaciers (0.01–0.5) dominate in terms of total number (749 glaciers), but the vast majority of the glacier area belongs to medium or 1.0–5.0 km² (354.0 km²) and large or > 10.0 km² glaciers (221.9 km²) (Table 2, Fig. 4b).

The mean elevation of the glaciers ranged from 3350 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (southern slope) to 3520 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (northern slope), with an average of 3475 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 5). Most of the glacier number (1476) and area (697 km²) in 2020 belong to north-facing glaciers (mean aspects N, NW, and NE), while relative area and number of E- and W-exposed glaciers are very small (Fig. 6).

Figure 7 (a) Glacier area vs. maximum and minimum elevation in 2020. (b) Glacier aspect vs. mean elevation in 2020. (c) Spatial distribution of mean elevation (colour-coded) for all glaciers in the Greater Caucasus larger than 0.1 km² in 2020. (d) Close view of the central part of the Greater Caucasus.

The glacier termini are located around an average minimum elevation of 3159 m a.s.l while the average maximum elevation is 3561 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Consequently, large valley glaciers have lower termini, while smaller glaciers have higher snout positions. All other topographic parameters (e.g. maximum, minimum, and mean elevations) depend on morphological type, aspect, and size class of the individual glaciers. Figure 7a and b show the glacier area distribution according to the maximum and minimum elevation and glacier aspect vs. mean elevation, while the colour-coded map at Fig. 7c shows spatial distribution of mean elevation for glaciers larger than 0.1 km² in 2020.

In 2020, the Elbrus Massif had a total glacier area of 107.7 ± 1.6 km² (Fig. S1). The three glaciers Bezengi (34.8 ± 0.8 km²), Karaugom (23.6 ± 0.3 km²), and Dzhikiugankez (19.4 ± 0.2 km²) are now the largest glaciers of the Greater Caucasus and are all located in Russia. Overall, there are 14 glaciers > 10 km² in the Greater Caucasus with a total area of 221.9 km². Three glaciers are situated in Georgia and 11 in Russia.

Figure 8 Greater Caucasus glacier area decrease by section and slope in 2000–2020.

5.3 Glacier change in 2000–2020

Results from our study on glacier area change indicate a significant decrease in the glaciers in the Greater Caucasus between 2000 and 2020 (Fig. 8). The total ice area loss between these two periods was 320.6 ± 45.9 km² or 23.2 ± 3.8 % (−1.16 % yr⁻¹). The eastern part experienced the highest absolute decrease of −1.82 % yr⁻¹, while the Elbrus Massif experienced the lowest rate of −0.57 % yr⁻¹. Compared to other sub-regions, the western region had also somewhat higher change rates (−1.45 % yr⁻¹). The Elbrus Massif has the largest glacier mean area, changing from 6.07 km² in 2000 to 3.98 km² in 2020.

Figure 9 Averaged annual area change rate (% yr⁻¹) from 2000 to 2020 for the seven glacier size classes in all sections and slopes of the Greater Caucasus.

The smallest size classes of glaciers (0.01 to 0.1 km²) experienced the highest area loss rates across all regions, with maximum rates in the eastern Greater Caucasus (Table 2, Fig. 9). The 0.1–0.5 km² size class also experienced high area loss rates (up to −2.9 % yr⁻¹ in the eastern part). For the larger size classes (> 1.0 km²) the loss rates are smaller and more similar. The difference in the loss rate between northern and southern slopes is not significant. Overall and similar to most other regions in the world, the observed relative area loss rates decrease towards larger glaciers.

Figure 10 Comparison of cumulative curves of terminus changes in 1960–2020: (a) for glaciers with size class > 10 km² and (b) for glaciers with size class 5–10 and 1–5 km². In both panels the dotted lines only connect the four measurement points. Data for 1960 and 1986 were taken from Tielidze and Wheate (2018).

From 16 selected glaciers (> 1 km²), the Lekhziri Glacier (43^∘9^′26^′′ N, 42^∘45^′54^′′ E) experienced the highest absolute retreat (1395 m or 69.8 m yr⁻¹) between 2000 and 2020, when the annual retreat of Lekhziri Glacier was ∼ 33 m in 1960–1986 and ∼ 13 m in 1986–2000 (Fig. 10a, Table S3 in the Supplement). Relatively small glaciers (1–5 and 5–10 km²) also experienced higher terminus retreat over the last 20 years, compared to previous time periods (Fig. 10b, Table S3). The smallest retreat between 2000 and 2020 from the selected glaciers was observed for Dolra Glacier (43^∘10^′10^′′ N, 42^∘31^′29^′′ E) with 178 m or 8.9 m yr⁻¹ (Table S3).

Table 3 Glacier area and count changes in the Greater Caucasus according to different inventories from 1911 to 2020.

Figure 11 Comparison of glacier area decrease rates in the Greater Caucasus for six different periods.

Figure 12 An example of a century-long area changes of Tsaneri Glacier (43^∘3^′25.68^′′ N 42^∘59^′1.92^′′ E) in 1890 (topographical map – X14, 1:42 000) (a); 1960 (topographical map – k_38_27, 1:50 000) (b); 1986 (Landsat 5 TM – 6 August 1986) (c); 2000 (Landsat 7 ETM+ – 5 September 2000) (d); 2014 (Landsat 8 OLI – 3 August 2014) (e); and 2020 (Sentinel 2B – 11 September 2020) (f). Outlines for 1890–2014 were taken from Tielidze (2016). Image credits: Panels (c), (d), (e): https://earthexplorer.usgs.gov (last access: November 2020); (f) Copernicus Sentinel data 2020.

6 Discussion

6.1 Comparison with previous investigations

In comparison to previous studies, our analysis reveals that the overall decline in glacier extent between 2000 and 2020 in the Greater Caucasus is 4 times higher than it was between 1911 and 1960, 3 times higher than it was between 1960 and 1986, and 2 times as high as it was from 1986 to 2000. An unprecedentedly higher decline was observed over the last 6 years, between 2014 and 2020 (Table 3; Figs. 11 and 12). Hence, our century-long comparison showed a clear decrease in glacier area in the entire region, which became much more pronounced over the last 20 years.

The observed glacier shrinkage in the Greater Caucasus from 2000 to 2020 (−1.16 % yr⁻¹) is similar to in the European Alps where Paul et al. (2020) reported a −15 % (or −1.3 % yr⁻¹) area reduction between 2003 and 2015/16. Direct comparisons with other glacierized regions are difficult because they are subject to different dynamics and size class distributions. Most of the related studies also cover different time periods. However, annual area loss rates larger than −1 % yr⁻¹ over the past decades have been reported for several regions in the world (e.g. Liu et al., 2020; Miles et al., 2021).

Figure 13 Comparison of glacier outlines from the RGI (v6) and GLIMS (in green) with the outlines from the new Caucasus glacier inventory (in yellow). (a) RGI nominal glaciers (circles) and glacier outlines derived during this study. The 28 July 2000 Landsat 7 image (Table S1) is used as the background. (b) The GLIMS outlines (an example of inconsistent registration) and glacier outlines derived during this study. The 12 August 2000 Landsat 5 image (Table S1) is used as the background. Image credit: https://earthexplorer.usgs.gov (last access: November 2020).

In comparison to existing glacier inventories, we found regionally large discrepancies that have now been corrected. The outlines included in the RGI v6 and GLIMS (2000) database were mostly created based on Landsat and ASTER imagery from 1999–2004 (Pfeffer et al., 2014; Khromova et al., 2009; Khromova et al., 2016). By detailed visual inspection, we found partly large differences between RGI v6, GLIMS (2000) outlines, and our database that was compiled using Landsat imagery from 2000. The RGI v6 contains nominal glaciers (circles) in the eastern and western Greater Caucasus, as well as the side ranges in the central Greater Caucasus that were replaced by real glacier outlines in our study (Fig. 13a). The GLIMS outlines also involve a horizontal geolocation shift (Fig. 13b), which appears to be associated with a shift in the ASTER images used (Tielidze and Wheate, 2018).

Figure 14 Comparison of cumulative glacier areas sorted for seven size classes for the RGI (v6), GLIMS, and current glacier inventory from 2000.

The RGI v6 contains 1638 glacier outlines with a total area of 1276.9 km². This is 548 fewer glaciers and ∼ 105 km² (∼ 7.5 %) less glacier area than mapped for this inventory. The largest differences were found for glaciers in the size class 1–5 km² (Fig. 14). The GLIMS database for the Caucasus region contains an even smaller number and area of glaciers than in the RGI (v6). In particular GLIMS does not contain the majority of glacier outlines from the eastern Greater Caucasus, resulting in 891 fewer glaciers and ∼ 270 km² (∼ 19.5 %) less glacier area than in our new database.

6.2 Uncertainties and limitations

The uncertainty of the mapping was assessed by a comparison of glaciers derived from multiple digitizations by different operators and using the buffer method. The resulting average uncertainty was less than ∼ 5 % of the mapped area, confirming the uncertainty estimate for the entire Greater Caucasus based on the buffer method (∼ 4 %). The major sources of uncertainty include the correct interpretation of debris cover, seasonal snow, and shadows, which can all impede accurate glacier mapping. Using imagery from a different date and local knowledge, the debris cover and shadow error have been partly resolved for some glaciers; while incorrect identification of seasonal snow generally affects small glaciers more than larger ones, where possibly included snow fields do not make up a large percentage of the total area.

We have not analysed here the temporal evolution of debris-covered glacier parts, as this is considerable extra effort and because we wanted to keep this study focused on the new inventories and the change assessment. However, we intend to analyse changes in debris cover in a separate study that might also consider recently developed methods (Holobâcă et al., 2021).

We used the ASTER GDEMv3 from around 2010 to derive topographic information for each glacier, although it fits neither 2000 nor 2020. The simple reasons for this decision are larger artefacts found in the SRTM DEM from 2000 and that a year 2020 DEM was not available for the study region. Accordingly, for shrinking and retreating glaciers mean and median elevations are underestimated for 2000 and – along with minimum elevation – overestimated for 2020. We assume that the related biases are within the uncertainty of the GDEM for most glaciers, but we wanted to stress that they have to be considered when working with the data. Unfortunately, this caveat is common in most similar studies (e.g. Paul et al., 2020) as the repeat frequency of freely available DEMs is still small. The impact of the wrong DEM timing on mean slope and aspect should be negligible.

Figure 15 Changes of mass balance of (a) Djankuat and (b) Garabashi glaciers in 1967–2019 and 1983–2019, respectively (WGMS, 2022).

6.3 Climatic and mass balance trends

Temperature data from Terskol meteorological station (northern Greater Caucasus – 43^∘15^′29^′′ N, 42^∘30^′51^′′ E) (Fig. 1) indicate annual air temperature increase by ∼ 1 ^∘C (from 11.5 ^∘C to 12.5 ^∘C) during the summer period (June, July, August) in 2000–2019 in contrast to a decreasing trend in October–May precipitation at the same time (from ∼ 720 mm to ∼ 650 mm) (Rototaeva et al., 2019). Increased summer temperature was also observed at the Mestia meteorological station (southern Greater Caucasus – 43^∘2^′56^′′ N, 42^∘44^′17^′′ E) (Fig. 1) between 2000 and 2014 (Tielidze et al., 2020c). Furthermore, the extension of ablation season over the last 2 decades was confirmed by instrumental measurement from both the northern (Garabashi Glacier) and southern (Zopkhito and Chalaati glaciers) Greater Caucasus (Rototaeva et al., 2019; Tielidze et al., 2020c).

The increased temperatures are also reflected in mass balance observations of two WGMS reference glaciers in the Caucasus region such as Djakuat (43^∘11^′48^′′ N, 42^∘45^′28^′′ E) and Garabashi (43^∘18^′15^′′ N, 42^∘28^′5^′′ E). They both show strong negative mass balances between 2005 and 2019 (Kutuzov et al., 2019; Rets et al., 2019; WGMS, 2022), resulting in a much higher ice loss in this time period than accumulated before 2005 (Fig. 15). Furthermore, assessment of glacier mass changes in the Caucasus region using the geodetic method over the period 2000–2019 (Hugonnet et al., 2021; Tielidze et al., 2022) showed a 3-fold increase in the rate of glacier mass loss.

It might be possible that the increase in incoming short-wave solar radiation in the high Caucasus mountains observed since the 1980s (10 W m² over 10 years) has played a significant role in the accelerated mass loss of glaciers in recent years. 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 (Toropov et al., 2019). Moreover, a decrease in the albedo of the glacier surfaces due to an increase in the concentration of mineral particles can be another possible reason of amplified glacier mass loss. Two different pollution events (5 May 2009 and 23 March 2018) are especially noteworthy, when an extreme amount of dust from the Sahara was deposited on the Caucasus glaciers, which sharply changed the albedo and accelerated melting in the accumulation areas (Kutuzov et al., 2013; Dumont et al., 2020). Due to additional factors involved for area changes (response times, ice thickness distribution), we do not relate here the observed more negative mass balances with the increased area loss over the same time period. However, in particular for the thin ice near the glacier terminus, we cannot exclude that the strong recent mass loss also contributed to the increased area loss.

7 Conclusions

We have presented the new Caucasus glacier inventory derived from manual delineation of glacier outlines based on medium-resolution (Landsat, Sentinel) and high-resolution (SPOT) satellite imagery acquired around 2000 and 2020. Within the entire Greater Caucasus, the total glacierized area mapped for 2000 and 2020 is 1381.5 ± 58.2 and 1060.9 ± 33.6 km², respectively, resulting in 23.2 ± 3.8 % (320.6 ± 45.9 km²) or −1.16 % yr⁻¹ reduction in total glacier area over the last 20 years. Glaciers < 0.5 km² contributed nearly 35 % to the total area loss but covered only 17 % of the total area (in 2000).

Glaciers in the western Greater Caucasus mostly have a lower mean elevation compared to glaciers in the central and eastern sections, indicating decreasing precipitation amounts from west to east. The highest area loss was observed in the eastern section, which is likely related to the decreasing glacier size to the east as relative area change rates increase towards smaller glaciers. The lowest decrease rate in the entire region was observed on the Elbrus Massif, which can be explained by the largest glacier area class dominating and maybe also an elevation that is sufficiently high to accumulate solid precipitation.

A century-long comparison with glacier areas mapped in previous inventories reveals a strong increase in area loss rates, nearly 4 times higher in 2000–2020 than it was between 1911 and 1960. Combined with the recent dominance of strongly negative mass balances, it can be expected that the Caucasus glaciers will continue to decline in the future under current climatic conditions.

With these two new glacier inventories for the Greater Caucasus, we have corrected errors in previously available datasets and hope that they will improve our understanding of climate change impacts at a regional scale and support related modelling studies by providing high-quality validation data.