Recent retreat of the Elbrus glacier system

Abstract

The glacier system covering Europe's highest mountain, Elbrus, has exhibited an accelerated retreat since 1980. Some studies have related this retreat to a significant summer temperature increasing trend. Relief- and aspect-related parameters for the glacierized area have an important impact on glacier changes. In this paper, the changes in glacier area are identified, quantified and correlated with relief parameters for the period 1985–2007. Spatial analysis was performed using the GLAM-CD (Glacier Mapper – Change Detector) algorithm. The input data for this algorithm were Landsat 5 images, the Aster Global Digital Elevation and the glacier outlines from the GLIMS project (Global Land Ice Measurements from Space). Regression analyses between glacier area losses and relief-related parameters indicate a significant positive relation with the altitude and a significant negative relation with the glacier surface area. In this context, we used a correlated component regression to model these relations. The model explains >50% of the total variation.

1. INTRODUCTION

The decline of mountain glaciers since the Little Ice Age (LIA) is an important topic in glaciology and related sciences. Mountain glaciers are vital sources of water for many regions in the world. Also, the risks associated with rapid ice melting can affect local communities.

Because glacier mass balance is very sensitive to climate variability, the dynamics of glaciers is a key indicator of low-frequency climate change (Oerlemans, 2005; Chinn and others, 2005). Worldwide, glaciers have retreated during the 20th century in response to global warming. Many regional analyses indicate a strong specific negative mass balance since 1990, with rapid losses especially in Patagonia, Alaska (Dyurgerov and Meier, 2005) and the European Alps (Zemp and others, 2005).

Caucasus glaciers have retreated since the mid-19th century (Solomina, 2000). Zolotarev and Kharkovets (2010) analysed the evolution of Elbrus glacierization since the LIA using cartographic and proxy methods, providing valuable information about past retreat rates.

The beginning of the 1980s was a landmark in glacial evolution in the Caucasus. Regional studies reported an accelerated retreat of Greater Caucasus glaciers since 1980 (Popovnin and Rozova, 2002; Stokes and others, 2006). Some studies related this retreat to an increase in summer temperature (Shahgedanova and others, 2005) or to a significant summer warming trend and a decline in winter precipitation (Holobâcă, 2013). Other studies highlighted potential consequences of the rapid glacial melt, such as glacial surges (Kotlyakov and others, 2002), multi-phase mass movements (Petrakov and others, 2008), and catastrophic detachment and high-velocity long-runout flow (Evans and others, 2009).

In a previous paper (Holobâcă, 2013), we found that climatic variability imposes the general trend for the Elbrus recent glaciation but relief parameters such as altitude, slope and aspect are very important for the glacial dynamics at the local scale. In this paper, we make a further contribution by relating the identified changes to relief parameters. Another important aim of this paper is to model these relations for future applications.

1.2. Study area

Situated in the northwestern part of the Greater Caucasus Mountain, Mount Elbrus is as an excellent candidate for local scale analyses. A well-developed ice cap covers the top of a dormant stratovolcano with two cones (Fig. 1). Mount Elbrus (5642 m), the highest peak in the Caucasus and Russian Federation, is covered by ~110 km² (in 2007) of ice and perennial snow (Table 1). The ice thickness surpasses 240 m on the Western Plateau (Mikhalenko and others, 2008).

Fig. 1. Mount Elbrus glacier system outlines, superposed over the 2007 colour composite of Landsat bands 5, 4, 3.

Table 1. Change of the glaciers surface area, relative losses and rate of change (1985–2007)

Notes: The average values for small glaciers (surface area <1.5 km²), medium glaciers (1.5 km²< surface area>10 km²) and large glaciers (surface area >10 km²) are printed in italic bold; the average values for Elbrus area are printed in bold.

According to the World Glacier Monitoring Service (WGMS) inventory, there are 24 valley glaciers radiating from this ice cap. The surface area of valley glaciers ranges from <0.1 to >22 km² (Table 1). The mean altitude of the Elbrus glacier system is about 4000 m, but the valley glaciers descend down to 2700 m in the southern and southeastern part (Table 2).

Table 2. Change of glacier altitude (1985–2007)

Notes: Median values for small glaciers (surface area <1.5 km²), medium glaciers (1.5 km²< surface area>10 km²) and large glaciers (surface area >10 km²) are printed in italic bold; the average values for Elbrus area are printed in bold.

The Elbrus valley glaciers, like other Caucasus glaciers, are affected by surge episodes (Kotlyakov and others, 2002). The associated rapid increases in ice flow rates can generate catastrophic detachments (Evans and others, 2009), and can be a major risk for downstream human structures and habitat. In addition, the Elbrus glacier system is an important source of freshwater for the North Caucasus region. Ice melt feeds the headwater of two major rivers, Kuban and Terek.

2. METHOD AND DATA

2.1. Spatial data

One of the major issues for our multitemporal surface analysis was to ensure the comparability of data. Only representative spatial and temporal data have been chosen. The length of the study period (1985–2007) was also an important concern. The recession trend observed for all Elbrus glaciers during the analysed period is consistent with climatic warming trends. The study period of 22 years appears to be long enough to avoid complications caused by non-climate related glacier perturbations (e.g. surge episodes). Also, Zolotarev and Kharkovets (2010) used a time period of 22 years (1957–79) to study Elbrus glaciers. This offers a good opportunity for comparison of results.

The Landsat 5 satellite images used in this study were acquired by the Thematic Mapper (TM) sensor, with 30 m equivalent pixel size. Both images were obtained from the end of the ablation season (late-July 1985 and early August 2007), when the maximum area of snow-free ice is exposed (Fountain and others, 1997). The scenes have minimal snow and cloud cover and were acquired before the first snowfall of winter. The source of satellite data is the United States Geological Survey Global Visualization Viewer (GLOVIS) Landsat archive (http://glovis.usgs.gov/). The orthorectified satellite scenes were co-registered using the ERDAS 9.1 software (Intergraph Corporation, Madison, AL, USA). The horizontal RMS error was <3.0 m.

Ice divides were located using Global Land Ice Measurement from Space (GLIMS) project data. This database is accessible at http://glims.colorado.edu/glacierdata/, and stores multitemporal glacier outlines and metadata (processing technics, scalar data and related literature) (Raup and others, 2007).

The relief parameters (altitude, slope and aspect) used were derived from an Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global DEM. The main advantage of this DEM, available at http://www.ersdac.or.jp/GDEM/E/ is that it has the same 30 m spatial resolution as the Landsat 5 TM images. The DEM was registered using heights from 1:50 000 topographic maps, accessible on the World Wide Web at http://mapstor.com/mapsets/russia-maps.html and 22 ground control points (GCPs) obtained with a GPS receiver in the field. The registration of the DEM was made using ERDAS 9.1 software (Intergraph Corporation, Madison, AL, USA). The method uses a second-order polynomial transformation, which allows a shift across the image that follows the curve of the second-order polynomial. The horizontal RMS error was in the range of 15 m (0.5 of DEM resolution).

2.2. Spatial analysis

Spatial analysis was carried out in multiple steps (Fig. 2): (1) derivation of the Elbrus glacier outlines from multitemporal satellite scenes using the Glacier Mapper – Change Detector (GLAM-CD) algorithm; (2) definition of the ice divides for individual glaciers from GLIMS data; (3) identification of the glacier outlines for individual glaciers combining two categories of spatial data (glacier–non-glacier outlines and ice divides); and (4) identification of the spatial patterns of the glacier changes.

Fig. 2. Automatic extraction and surface analysis algorithm.

2.2.1. Derivation of the Elbrus glacier outlines

The GLAM-CD is a logical structured algorithm that allows evaluation of relief-related changes in the glacial area from two time-comparable satellite images (Holobâcă, 2013). This ARC GIS toolbox uses an image ratio method in order to enhance the variation in the slope of the spectral reflectance curves between two different spectral ranges. The GLAM-CD has two tools, one for proper ratio images and the other for the Normalized Difference Snow Index (NDSI) calculation. The NDSI, a widely used index for automated glacier delineation (Paul and others, 2002), was employed in this analysis. As snow strongly reflects visible light (e.g. green) but absorbs middle infrared light, NDSI is high over snow and ice. Acceptable snow cover results can be found with NDSI thresholds in the range, 0.25–0.45 (Hall and others, 1995). In this study, glacier outlines are derived using a threshold value of 0.4, which provides better results for glacier delineation.

We use a semiautomatic method to extract glacier outlines from the Landsat 5 images. The initial glacier outlines obtained with the GLAM-CD tool were checked visually using the colour composite of Landsat bands 5, 4, 3 and 3, 2, 1. This visual correction for classification errors, debris covered ice and shadow analysis is important due to the limitations of full automated methods (Paul and others, 2013).

2.2.2. Definition of ice divides for individual glaciers

Following algorithms developed by Manley, GLIMS uses watershed commands in ArcInfo to separate individual glaciers from large contiguous ice masses (Raup and others, 2007). In this multitemporal study, the same ice divides have been used for individual glacier definition. IDs or names from the GLIMS database were adopted for glaciers rather than the local name from topographic maps. In order to avoid misidentifications we also employed the WGMS IDs. The ice divides for individual glaciers were checked visually, and manually improved using colour composite Landsat bands (5, 4, 3 and 3, 2, 1). Special attention was given to internal rocks in the ice divide area.

2.2.3. Identification of glacier outlines for individual glaciers

The individual glacier outlines have been obtained combining two categories of vector spatial data. These vector data were the Elbrus glacier outlines and the individual glacier ice divides.

2.2.4. Identification of the spatial patterns of glacier changes

We analysed the vector files from the two-time positions of the glacier outline to identify changes in the glacierized area.

The relative losses of individual glaciers were correlated with the characteristics of the relief parameters (Table 1). We then tested the correlation coefficient at the 0.05 significance level. In addition, a Correlated Component Regression Linear Model (CCR.LM) was employed to model the relation between glacier changes and relief parameters. This regression method uses fast cross-validation to determine the amount of regularization in order to produce reliable predictions from data with correlated explanatory variables. We use XLSTAT software (Addinsoft, USA) to build the model.

3. RESULTS

Recession of glaciers was the main characteristic during the analysed period. The size of the change varied dramatically from glacier to glacier (Tables 1 and 2). Climate is thought to be the main driver for the changes, but local conditions can alter the response of individual glaciers to the warming trend. In order to better understand these complex relations, a quantitative approach is presented below.

3.1. Glacier areal change

Table 1 indicates that glacier size is important for the intensity of the recession. The statistical relation is negative and significant at 0.05 p-level (Table 3). Small glaciers, with an area <1.5 km², dramatically diminished when considering relative changes. The relative losses for these glaciers with an average value of 32.2% and ranging from 12 to 64%, surpassed the average value for the Elbrus glacier system as a whole (8%).

Table 3. Correlation coefficient values (R) between glacier losses and relief parameters

Notes: Correlation values statistically significant at 0.05 level of signification are printed in bold, min, minimum; max, maximum; S25%, surface slope averaged over the lowermost 25% of the glacier; E, easting; N, northing.

For the most part, for the medium-sized glaciers, with an area ranging from 1.5 to 10 km², glacier losses were <8%. The average value for the medium glaciers was 7.8%. The most frequent range of values was 2.5–3.6%. We noticed two major exceptions: Garabashi Glacier, with the smallest relative loss value within the Elbrus glacier system (0.8%) and Glacier 576 (SU4H08004316), which lost 39% of its surface area in the analysed period. The Garabashi Glacier is situated in the southeast part of Mount Elbrus. The redistribution of snow from the northwest to the southeast side of the mountain range by the dominant north-westerly winds can explain this small glacier relative loss.

The values of relative losses for glaciers with an area >10 km² were greater than expected with an average value of 7.0%. Dzhikiugankez Glacier, the most extensive glacier in the study area (22 km²), lost 19% of its surface area. The relative loss values for Ulluchiran (5.7%) and Bol'shoy Azau (9.3%) glaciers were close to the value for the overall glacierized area. Only Malyi Azau Glacier had a relative loss in the expected range (1.8%). The important relative losses for the Dzhikiugankez Glacier can be explained by the role of the intense post-volcanic activity, especially the influence of thermal and fluid flows in the north-eastern part of the Elbrus volcano (Masurenkov and Sobisevich, 2012).

3.2. Glacierized altitude change

Many glaciological studies suggest that altitude is the most important independent variable and is a fundamental control of both temperature and precipitation (Evans, 1969; Sugden and John, 1969; Flint, 1971; Price, 1981). Changes in altitude also indicated a generalized glacial retreat for all the Elbrus glaciers. The results had some spatial patterns, exhibited in Table 2.

Change in minimum altitude is an objective indicator of frontal glacier retreat. The largest changes, ranging from −93 to −131 m, characterized the glaciers that descended below 3000 m. The loss in minimum altitude for the glaciers situated above 3000 m was complex. The loss values greater than the Elbrus average change (−47 m) were observed for medium and large glaciers that descended at low altitude (below 3200 m). The average values for these glaciers were −46.1 and −71.4 m, respectively. For small glaciers, altitudinal losses >47 m have been observed over a large altitudinal range (3532–3743 m). The medium range of values (−43 to −49 m) characterized the small glaciers with an average value of −34.5 m. The exception was the medium Karachaul glacier, which lost 45 m in minimum altitude. Minor losses in minimum altitude, sometimes below the DEM resolution, characterized the smallest glaciers, but were also observed for some medium and large glaciers (Table 2).

The average altitude of the Elbrus glaciers rose during the analysis period. The average altitudinal change for the Elbrus glacier system was −29 m. Changes >−30 m characterized the glaciers with large area losses (Table 1) having west, northwest, southwest and northeast-facing slopes (Fig. 1). Bityuktyube Glacier had the biggest change in average altitude (−81.9 m).

For the most part, the medium range of average altitude change (from −20 to −30 m) characterized the medium and large glaciers with minor glacial area losses. The only exception was the 310 glacier with an area of 0.1 km².

The class of minor average altitudinal change values ranged between −0.4 and −20 m. Only three glaciers (532v, Irik and Garabashi) with an area >1 km² had minor average altitudinal losses. All these glaciers are situated in the south-eastern part of the Elbrus massive. On the other hand, the glaciers with area <1 km² had the smallest average altitude losses.

Fewer changes were observed at maximum altitude. Most glaciers had no change in maximum altitude during the analysis period (Table 2). Some of them experienced a decrease in altitude. The maximum decrease (−69 m) was observed for the 565a glacier. Small increases in maximum altitude, ranging from 2 to 17 m, were identified for three glaciers. The changes in the maximum altitude seem to be random and uncorrelated with other relief parameters.

3.3. Modelled relations between relief parameters and glacier losses

Table 3 presents correlation coefficients between glacier loss and relief parameters. The relative glacier losses were correlated with initial relief parameters for individual glaciers, e.g. surface altitude, slope and aspect (easting and northing). Correlation coefficients were tested at the 0.05 significance level.

As expected, glacier surface area has a negative correlation with glacier loss. Significant negative correlations between these two parameters were found during the studied period (Table 3). Generally, the small glaciers have been more affected by glacial retreat. Also, we found significant correlations between glacier loss and altitude parameter. Our analysis revealed a significant positive correlation between glacier loss and minimum altitude. This relation can be explained by the position of the small glaciers (with high loss) situated at higher altitude than the medium or larger glaciers. The medium or larger glaciers, particularly the glaciers situated in the south or south-eastern part of the Elbrus massive had lower minimum altitudes and minor glacier losses (Tables 1 and 2).

Change in glacial surface area and the average and maximum altitude are negatively correlated. The glacier losses were higher at low altitudes and the glaciers with high average or maximum altitude have been less affected by deglaciation. Correlation coefficients, ranging from −0.5 to −0.6, were statistically significant.

The difference between maximum and minimum altitude has a significant negative relation with relative glacier loss. This variable had the highest correlation coefficient value (−0.7) indicating a higher climatic sensitivity of the glaciers with smallest altitudinal range.

Of the parameters tested against glacier loss, the lowest (magnitude) correlation coefficients, ranging from −0.3 to 0.2, were observed for maximum slope and average slope. Regression of the surface slope averaged over the lowermost 25% of the glaciers indicated a correlation coefficient, 0.2 which was not statistically significant. On the contrary, the minimum slope had a significant correlation with glacier loss. In our opinion, this statistically significant correlation can be explained by the size of the glaciers. The medium and large glaciers with minor glacier losses had minimum slope values of zero (or very close to zero) and the small glaciers with significant glacial surface area losses had larger minimum slope values.

The aspect parameters included in the model were easting and northing. The easting had a direct relation with relative glacier losses (0.3), but the relation between changes in glacial surface area and northing becomes negative (−0.3). The results seem to indicate a preferred sheltered orientation of the glaciers from direct solar radiation.

A correlated component multiple linear regression model was built in order to relate the relief parameters to glacier losses. The model uses 11 variables, listed in Table 3. The dependent variable is relative glacial loss for individual glaciers.

The model equation is:

(1) $$\eqalign{{\rm PGL}(\% ) & = 113.34 - 10^{ - 6 \ast} S1985 + 0.005^{\ast} {\rm altmin} \cr &\quad - 0.004^{\ast}{\rm altmax} - 0.027^{\ast}{\rm altaver} \cr &\quad + 0.0006^{\ast}({\rm altmax} - {\rm altmin}) + 5.59^{\ast}{\rm slopemin} \cr &\quad + 0.11^{\ast}{\rm slopemax} + 0.49^{\ast}{\rm slopemed} \cr &\quad - 0.50^{\ast}S25\% + 6.34^{\ast}N - 0.42^{\ast}E} $$

where PGL is relative glacier loss, S1985 is the surface area (km²), altmin is the minimum altitude (m), altmax is the maximum altitude (m), altaver is the average altitude (m), slopemin is the minimum slope (°), slopemax is the maximum slope (°), slopemed is the average slope (°), S25% is the surface slope (°) averaged over the lowermost 25% of the glacier, N is northing and E is easting.

The standardized (β) and unstandardized (B) coefficients of the multiple regression are presented in Table 4. The standardized coefficients indicate which of the independent variables has greater effect on the dependent variable in a multiple regression (Schroeder and others, 1986). The results indicate that average altitude (β = −0.39) and minimum slope (β = 0.43) have the major impact on the relative glacier losses.

Table 4. Unstandardized (B) and standardized (β) coefficients of the multiple regression

Notes: min, minimum; max, maximum; S25%, surface slope averaged over the lowermost 25% of the glacier; E, easting; N, northing.

The regression does a good job of modelling glacier losses. The adjusted r ² has a value of 0.67. Therefore, more than half the variation in relative loss is explained by the model. Also, the cross-validation value of the model is 0.14. Two usual tests for heteroscedasticity (Breusch and Pagan, 1979; White, 1980) were performed on model residuals. The significance level was 0.05. The results of both tests indicate that the null hypothesis (H ₀: Residuals are homoscedastic) cannot be rejected at the selected significance level.

4. DISCUSSION

Our technique, combining semiautomatic extraction of glacier outlines with surface analysis has some major advantages. Automatic delineation of glacier outlines on a TM scene is up to 300 times faster than manual techniques (Paul and others, 2013) and even our semiautomatic approach is much faster. In addition, there is apparently little loss of accuracy associated with automatic glacier outline delineation. Relative area differences between automatic and manual delineation techniques range from 2 to 5% (Paul and Kääb, 2005; Bolch and Kamp, 2006; Andreassen and others, 2008). In our case, differences fall within this range. Larger values have been observed for large and complex glaciers (e.g. Dzhikiugankez glacier – 3.4%, Irik glacier – 2.1%, Kyukyurtlyu glacier – 1.7%, etc.), while the difference was null or insignificant (<0.5%) for the small glaciers.

The integration of individual glacier ice divides with the GLIMS database considerably reduces processing time. Delineation of ice divides requires multiple spatial data sets (detailed DEM, high definition remote sensing data, aerial photography, field survey, etc.) and uses a complex algorithm. Use of the GLIMS ice divides also assures compatibility of multitemporal data. In fact, the reduction, or even elimination, of consistency problems is one of the most important GLIMS aims (Raup and others, 2007). Also, automatic standard algorithms can facilitate meaningful comparisons between multitemporal analyses within the same region (Bishop and others, 2004).

Debris-covered areas of glaciers remain the most important source of systematic error in automatic glacier outline delineation. This issue can be solved using semi-automated methods proposed by Paul and others (2004), utilizing complex processing techniques and an accurate DEM. The obtained data are manually edited. More recent, accurate debris-covered glacier data have been acquired using summer coherence images from microwave sensors (Atwood and others, 2010; Frey and others, 2012).

Spatial analyses have been designed to detect changes in glaciated areas and to correlate these changes with relief parameters. Our findings support previous studies (e.g. Popovnin and Rozova, 2002; Shahgedanova and others, 2005; Stokes and others, 2006 etc.) that reported accelerated retreat of glaciers in the Caucasus Mountains. This is more evident when the relative rate of change of glacier surface area (% a⁻¹) in the analysis period is compared with the Zolotarev and Kharkovets (2010) data (Fig. 3). The relative rate of change was calculated for three previous periods: 1850–87, 1887–1957 and 1957–79 and computations were made employing the absolute values of the glacier surface area change for the total Elbrus glaciation (source of data: Zolotarev and Kharkovets, 2010, p. 24, Table 3).

Fig. 3. Retreat trends in Elbrus Mountain since 1850 (% a⁻¹). Data source: Zolotarev and Kharkovets (2010); p. 24, Table 3

The relative rate of change of glacier surface area in our analysis period exhibits the largest decrease since the LIA in Elbrus Mountain. The rate of change of glacier's surface area from 1985 to 2007 for the Elbrus glacier system (−0.37% a⁻¹) is more than double that during the period 1957−97 (−0.16% a⁻¹). The actual rate is even greater when compared with the change rate of −0.15% a⁻¹ during the period 1887–1957 and almost double the rate of −0.20% a⁻¹ observed in the 1850–87 period (Fig. 3).

Surface analysis of glacial changes confirms the findings of previous studies (e.g. Shangguan and others, 2006; Lemke and others, 2007; Frauenfelder and Kääb, 2009; Bolch and others, 2010; Le Bris and others, 2011 etc.), that small glaciers (area <1.5 km²) and low glaciated areas (below 3200 m) are most affected by deglaciation. The major exception was Dzhikiugankez glacier, the larger Elbrus glacier, which exhibits high loss. This exception may, in part, be due to factors other than climate or relief, e.g. endogenic factors. Masurenkov and Sobisevich (2012) found that intense post volcanic activity, especially the influence of thermal and fluid flows, was responsible for glacial loss in the north-eastern part of Mount Elbrus.

The results of the regression analysis indicated statistically significant relationships between relative glacier losses and glacier surface area and altitude. Correlations with the average and maximum slope and with easting and northing of the aspect were not statistically significant. Only the minimum slope had a statistically significant relation with the melted glacial surface area. This situation may indicate the state of the dynamic equilibrium of the glaciers. The glaciers loosing area least have near zero minimum slope values, while the glaciers affected more by deglaciation have non-zero slope values.

Previous studies indicate some issues when using relative area changes as a parameter for inventory change (Nuth and others, 2013). The most important reported problems are the dependence of this parameter on the terrain geometry and heteroscedasticity of the variability for the glaciers having different sizes. In our approach, the correlated component regression linear model related the parameters of the relief with the relative area change. Also, the heteroscedasticity test revealed that the residuals of the model are homoscedastic.

5. CONCLUSION

The results of the change analysis indicate a general glacier recession trend of 0.4% a⁻¹ during the studied period. This recession is very uneven among the Elbrus glaciers. In our opinion, this situation is caused by the impact of local geographical factors superposed on the general climatic driver of change.

A technique that combines a semiautomatic extraction of glacier outlines with a surface analysis is a fast and consistent method of detecting glacial change. These qualities can be very useful if the analysis is extended to the regional scale.

We developed a regression model between relative glacier loss and the relief parameters, which explains more than 67% of the total variation. These results can help in water supply planning and management of rapid glacial melt-associated risks in the Mount Elbrus area. It can also lead towards further models that can predict future glacial area evolution at both local and regional scale.