Post-Little Ice Age Glacier Recession in the North-Chuya Ridge and Dynamics of the Bolshoi Maashei Glacier, Altai

Abstract

The glacier recession of the North-Chuya ridge, Altai, after the maximum of the Little Ice Age (LIA) is estimated based on remote sensing and in situ studies of the Bolshoi Maashei glacier. The glacier area decreased from 304.9 ± 23.49 km² at the LIA maximum to 140.24 ± 16.19 km² in 2000 and 120.02 ± 16.19 km² in 2021. The average equilibrium-line altitude (ELA) rise after the LIA was 207 m. The reduction of glaciers was caused by the warming trend, most rapid in the 1990s, and by the decrease in precipitation after the mid-1980s. The volume of glaciers decreased from approximately 16.5 km³ in the LIA maximum to 5.6–5.8 km³ by 2021. From the LIA maximum to 2022, the Bolshoi Maashei glacier decreased from 17.49 km² to 6.25 km², and the lower point rose from 2160 m to 2225 m. After the LIA, the glacial snout retreat was about 1 km. The fastest retreat of the glacier terminus was estimated in 2010–2022 as 14.0 m a⁻¹ on average. The glacier mass balance index was calculated, with the results showing a strong negative trend from the mid-1980s until now. Strong melt rates caused the increase in the area of the Maashei lake, which could lead to the weakening of its dam, and prepared for its failure in 2012. The current climatic tendencies are unfavorable for the glaciers.

1. Introduction

Glaciers are natural indicators of climate change of a different scale [1,2,3]. Their indicator role is higher for mountainous regions, where meteorological stations are scarce and smaller mountain glaciers respond more quickly to climate fluctuations due to their smaller size. The role of mountain glaciers is also great in a number of other processes of local and global scales: mountain glacier runoff is vital for the water supply of arid regions [4,5,6,7,8,9], its changes can cause the activation of glacial lakes’ outbursts and other dangerous exogenous processes [10,11,12,13], and the mountain glacier runoff is the main contributor to the changes of the sea level [14,15,16,17].

Climate change has a major impact on the mountain glaciers that can influence all of those processes. Measurements show that the Earth’s global average near-surface temperature has increased by about 0.8 °C since the 19th century [18]. This climate change has been accompanied by the worldwide general tendency of the glacier retreat after the Little Ice Age (LIA) maximum [19,20], that accelerated in the 21st century [16].

The current worldwide trend of retreat of the glaciers is also present in the Altai mountains, and this retreat has accelerated in the past several decades [21,22,23,24,25]. However, there is still little information about the shrinkage of the Altai glaciers after the LIA; on the contrary, most of it is about the largest glaciers, and this information is usually outdated due to the lack or reduction of monitoring of the retreat of the largest glaciers of Altai in the 1990s–2000s. Recently, we started working on this problem in order to fill the gaps in knowledge of the glacier dynamics in Altai after the LIA [26,27]. The aim of the present article is to reveal the LIA shrinkage after the LIA for the North-Chuya ridge, to reconstruct the retreat of one of the largest glaciers of the Altai–Bolshoi Maashei glacier (the systematic monitoring of its retreat has been stopped since 1980), and to identify the relationship between glacier dynamics and climate change.

2. Study Area

The Altai mountainous country is a territory located at the junction of Central Asia and the mountains of Siberia, which combines the features of their climate: the weakened influence of the oceans on the climate; the summer maximum precipitation; large annual temperature contrasts on the one hand; and western transport as the main source of precipitation, negative annual temperatures, and development of permafrost on the other. The territory of Altai is divided between four countries: Russia, Mongolia, China, and Kazakhstan. Russian Altai is the highest part of the Altai Mountain system, reaching 4506 m a.s.l. (Belukha Mountain). The combination of extensive planation surfaces with high-mountain alpine relief is a typical feature of the Altai relief. In its central part, mountain ranges and massifs rise to 3000–4000 m a.s.l. (Katunsky, South-Chuya, North-Chuya ridges, Figure 1). These ridges carry both modern glaciers of various morphological types, from large valley glaciers to small cirque and hanging ones, as well as traces of their larger extent in the past, including the LIA. Those traces are presented by such landforms as moraines, cirques, u-shaped valleys, and nival niches.

The North-Chuya ridge stretches about 100 km in the east-west direction. The ridge reaches its greatest height in the central part, known as the Bish-Iirdu mountain junction, in which the main glaciation of the ridge is concentrated. The average altitude of the ridge here is about 3600 m a.s.l., and a number of peaks exceed 4000 m a.s.l. (Maashei-bash: 4177.7 m a.s.l.; Aktru: 4044.4 m). This highest part of the ridge is the watershed of the Chuya river (its tributaries Maashei, Aktru, and Tete rivers on the north-eastern slope; Chagan-Uzun (Jelo river) on the south-eastern slope; and Argut river (south-western slope)) with its tributaries Shavly and Yungur rivers on the north-western slope, and Karagem on the south-eastern slope (Figure 2). All those rivers belong to the Katun river basin, one of the sources of the Ob river, belonging to the Arctic ocean basin.

Directly on the territory of the North-Chuya ridge there is only the Aktru weather station, which was founded in 1972, but its data are not available in the public domain, and in addition, observations at it were interrupted. The other three weather stations are located relatively close to the ridge. One of them, the Kosh-Agach weather station, is located in the center of the Chui basin. The use of data from this meteorological station is of little use in characterizing the climate of the North-Chuya range, since it is located in the orographic shadow of the North-Chuya and South-Chuya ranges in arid conditions (about 120 mm of precipitation per year), while the North-Chuya range, especially its northern slope, is relatively accessible to moisture-bearing western atmospheric flows. About 30–40 km to the west of the North-Chuya ridge, there are meteorological stations Akkem (2050 m a.s.l.) and Kara-Turek (2600 m a.s.l.). The high altitude of the weather stations, their proximity to each other, and their location in the direction from which the moisture-bearing flows come, suggest using their data for calculating temperature and precipitation gradients and to obtain a climatic characteristic of the glacial zone of the North-Chuya range. According to Kara-Turek data, which is the highest mountain weather station in the Altai region, the average multiannual (1961–2021) temperature is −5.6, the average multiannual summer temperatures (June, July, August) are 6.1 °C, and the average multiannual precipitation is 610 mm. In the glaciated zone, according to multiannual observations on the glaciers of the Aktru basin, annual precipitation is about 800–900 mm [29].

The glaciers of the North-Chuya ridge were first discovered in 1898 by Sapozhnikov who visited the North-Chuya range and discovered two glaciers in the Aktru river basin and the Dzhelo glacier [30]. In 1924, M.V. Tronov and B.V. Tronov explored the Bish-Iirdu mountain junction and obtained new information about the glaciers of the ridge. Moreover, in the river Maashei valley, they discovered a previously unknown glacier, called the Big Maashei Glacier of North-Chuya. The results were included in the first catalogue of the glaciers of Altai [31]: the first estimate of the total glaciation of the North-Chuya ridge therefore included 118 glaciers with an area of 127 km². The next investigation of the North-Chuya glaciers, which was made as a part of the work for the USSR Glacier Inventory, gave the following estimation: 168 glaciers with the total area of 143.16 km² [32]. Later, data from the Glacier Inventory were corrected and updated for 2003: 181 glaciers with a total area of 164.2 km² [33]. The latest estimate of the glaciation of the North-Chuya glaciers for 2017 gave an area of 112.9 km² [34].

Until 2012, in the valley of the Maashei river, at a distance of about 4.8 km, there was a lake with an area of approximately 0.25 square kilometers. The basis of the lake’s nutrition was glacial runoff. The dam of the lake was formed by a scree and a rock glacier. In 2012, as a result of heavy rains, the dam of the lake was destroyed, and the lake was drained catastrophically.

3. Materials and Methods

3.1. Field Glaciological and Glacio-Geomorphological Observations

We used field observation data to study the dynamics of the Bolshoi Maashei glacier. They include the results of observations by scientists who have visited the glacier since its discovery in 1924, and the results of our observations in September 2022.

Unfortunately, scientists did not directly observe the time of the LIA maximum in Altai. For the Katunsky glacier (about 45 km to the south west of the North-Chuya ridge), it is known that the glacier was advancing in 1835 [24,35], and by 1880 it had already retreated from its position in 1835 by 350–380 m [24]. In the North-Chuya range, observations started only after 1898 [30] when the glaciers were already retreating, and mass balance observations started in 1961 [36]. However, there were some years with a positive mass balance, especially between 1965 and 1995, and as a result the halting of the glacier retreat and even a slight advance of glaciers in the Aktru basin, most well studied by the moment, occurred in 1936, 1940, 1969, and 1993 (Right Aktru) and in 1911, 1936, 1960, 1966, 1979, and 1993 (Maliy Aktru) [37].

The Bolshoi Maashei glacier was first visited by B.V. and M.V. Tronovs in 1924, and then by M.V. Tronov in 1932 and 1937 [35]. In 1952, 1962, and 1975 the glacier was explored by P.A. Okishev [24,32], and further work was carried out by R.M. Mukhametov, who installed a benchmark near the edge of the glacier in 1980. Accordingly, for the period of 1924–1980, there are data on the retreat of the glacier edge based on direct measurements [24]. We have no way to evaluate the accuracy of these measurements (by the measuring tape) but based on our own experience of similar work and taking into account the shallow bottom sloping of the valley, we can assume that the error did not exceed 0.5 m and thus it was not taken into account in further calculations.

Since 1980, no systematic observations of the glacier’s dynamics have been carried out due to its difficult accessibility. At the same time, this glacier is one of the largest in Altai, and the originality of its dynamics was noted by M.V. Tronov, who indicated its stability and relatively slow reduction, which is associated with a vast firn basin lying at a high altitude, the bottom of which is much higher than the snow line (at an altitude of more than 3700 m with a snow line altitude of about 2900 m a.s.l.) [35].

In September 2022, we visited the Bolshoi Maashei glacier. During the survey, we performed GPS tracking of the glacial snout, found and geodesically surveyed the 1980 benchmarks of Mukhametov, and measured (both by measuring tape and by GPS) the distance from the benchmark to the glacial snout’s nearest point. We also created a new benchmark (27.8 m from the current position of the lowest point of the glacier), as the 1980 benchmark is about 340 m from the glacier, which complicates the measurement of retreat.

We also carried out geomorphological observations in the forefield of the glacier, which allowed the GPS to mark young, slightly vegetated moraines dating back to the period of glacier shrinkage after the Little Ice Age.

3.2. Remote Sensing

The remote sensing study was used for several purposes. The first group of tasks was to obtain information about the whole North-Chuya ridge (to delineate the present state of the glaciers, to reconstruct the positions of the glaciers in 2000, and to delineate the geomorphological features that are used to reconstruct the glaciers in the maximum of the LIA). The second group of tasks was to reconstruct the positions of the glaciers in the Maashei river basin for as many time points as possible in order to extend the glacier shrinkage data series.

We chose the images from the end of the ablation season, when the elevation of the snow line is maximal, and the images are with minimal snow cover (not after the snowfalls) and with low cloud cover. If clouds covered parts of the glacierized area in an individual image, another image from the nearest available date was collected.

For creating the present glacier inventory and reconstruction of the LIA glaciers, we used a Sentinel-2 image, and for the reconstruction of the glaciers in 2000, we used Landsat 7 imagery (

Table 1. Satellite imagery used in the study.

Date	Spacecraft	Spatial Resolution, m	Image ID
2021-09-08	Sentinel-2	10	S2B_MSIL1C_20210908T051649_N0301_R062_T45UWR_20210908T075210.SAFE
2019-08-02	Sentinel-2	10	S2A_MSIL1C_20190802T050701_N0208_R019_T45UWR_20190802T081536.SAFE
2016-07-01	GeoEye-1	1.65	10500100052CD500
2013-08-24	World View 02	1.8	1030010025D4D800
2010-09-20	LANDSAT_7	15.0	LE71440252010263ASN00 M
2008-08-29	LANDSAT_7	15.0	LE71440252008242SGS00 M
2007-08-27	LANDSAT_7	15.0	LE71440252007239PFS00 M
2006-09-09	LANDSAT_7	15.0	LE71440252006252PFS01 M
2004-08-18	LANDSAT_7	15.0	LE71440252004231ASN01 M
2000-08-07	LANDSAT_7	15.0	LE71440252000220SGS00 M
1993-08-19	LANDSAT_5	30.0	LT51450251993231ISP00 M
1993-08-12	LANDSAT_5	30.0	LT51440251993224ISP00 M
1989-08-17	LANDSAT_5	30.0	LM51440251989229ISP01 M
1980-07-03	LANDSAT_3	40.0	LM31550251980185AAA08 M
1980-08-08	LANDSAT_3	40.0	LM31550251980221AAA05 M
1968/08/10	CORONA	1.8	DS1104-1039DA004
1962/06/28	CORONA	1.8	DS009038052DF039

). For the reconstruction of the glaciers in 2000, we used a 2000-08-07 LANDSAT_7 image. Other images were used for the second group of tasks and also were used as an additional source of information for the first group of tasks.

Landsat images and CORONA images were provided by the USGS [38]. For Landsat, one common false-color composite (FCC) image combination was provided by TM bands 5, 4, 3, which help delineate clean glacier ice/snow and vegetation [25,39]. For Landsat 7 images, we used false-color composites (FCCs) that showed the differences in reflectance of landscape features. In particular, we used a Landsat ETM+ and TM 5, 4, 3 RGB composite (red: channel 5; green: channel 4; blue: channel 3), and the resulting effect was that snow and ice were clearly differentiated from clouds, debris, rock, or vegetation due to FCC image color differences [39]. Furthermore, Landsat 7 images taken by the ETM+ sensor included a panchromatic band 8, which was used for pan-sharpening to improve image resolution from 30 m to 15 m employing ESRI ArcGIS 10.4. We also used Landsat 2 MSS 765 and 654 RGB composites. We have developed a number of rules to avoid errors when identifying glacier boundaries (Table 2).

Table 2. Main mistakes in glacier delineation and the ways to avoid them.

Probable Mistake Description	Area and Length Error Sign	Ways to Avoid Mistakes
Debris-covered glacier edges misinterpreted as parts of the LIA complex or dead ice	-	Usage of active ice indicators: the “smooth” debris surface (Figure 3B), linear flow structures, constrained tributaries [40], and marginal water flows that join at the lowest point of the glacier [26]; in some cases, the debris-covered parts of the glaciers are not smooth, but are broken to blocks by crevasses (Figure 3A)
The areas of dead ice under the moraine cover delineated as parts of the glacier	+	Usage of dead ice indicators: surface subsidence forms related to downwasting of debris-covered dead ice [41], the rugged debris surface, melting ponds, unconstrained tributaries, pioneer plants [40], and water streams going into tunnels and going out of other tunnels down the slope [26] (Figure 3A)
Medial moraines on the surface of glaciers mistaken for rocky outcrops that separate glacial streams	-	Search for the areas where the moraine cover crumbles and the ice core of the moraine is exposed, which helps to correctly diagnose it
Seasonal snow cover areas, taken as parts of the glacier	+	Usage of additional imagery with less-extensive snow cover
Perennial snow patches taken for parts of the glaciers or small glaciers	+	Glaciers have an integral configuration, and snowfields often have openwork outlines in the plan; for glaciers, the images show ablation and accumulation zones; crevasses or bergschrunds visible on the surface of glaciers, as indicators of glacier movement
Shading of parts of the glacier and the adjacent non-glacial areas	-	Usages of additional images of the same area with different acquisition times and different angles of sunlight are compared; in other cases, there is adjusting the brightness and contrast of the image to make the shaded area visible, or usage of different channel combinations
Frozen glacial lakes taken for parts of a glacial tongue	+	Usage of additional images for the period when the ice on the lake has already molten; usage of digital elevation models, in which this area will look such as an absolutely flat territory, contrasting with the sloped areas of the glacial tongue
Icings	+	Icing ice is usually lighter than glacial ice due to its cleanliness; such areas appear flat on digital elevation models

Table 2. Main mistakes in glacier delineation and the ways to avoid them.

Probable Mistake Description	Area and Length Error Sign	Ways to Avoid Mistakes
Debris-covered glacier edges misinterpreted as parts of the LIA complex or dead ice	-	Usage of active ice indicators: the “smooth” debris surface (Figure 3B), linear flow structures, constrained tributaries [40], and marginal water flows that join at the lowest point of the glacier [26]; in some cases, the debris-covered parts of the glaciers are not smooth, but are broken to blocks by crevasses (Figure 3A)
The areas of dead ice under the moraine cover delineated as parts of the glacier	+	Usage of dead ice indicators: surface subsidence forms related to downwasting of debris-covered dead ice [41], the rugged debris surface, melting ponds, unconstrained tributaries, pioneer plants [40], and water streams going into tunnels and going out of other tunnels down the slope [26] (Figure 3A)
Medial moraines on the surface of glaciers mistaken for rocky outcrops that separate glacial streams	-	Search for the areas where the moraine cover crumbles and the ice core of the moraine is exposed, which helps to correctly diagnose it
Seasonal snow cover areas, taken as parts of the glacier	+	Usage of additional imagery with less-extensive snow cover
Perennial snow patches taken for parts of the glaciers or small glaciers	+	Glaciers have an integral configuration, and snowfields often have openwork outlines in the plan; for glaciers, the images show ablation and accumulation zones; crevasses or bergschrunds visible on the surface of glaciers, as indicators of glacier movement
Shading of parts of the glacier and the adjacent non-glacial areas	-	Usages of additional images of the same area with different acquisition times and different angles of sunlight are compared; in other cases, there is adjusting the brightness and contrast of the image to make the shaded area visible, or usage of different channel combinations
Frozen glacial lakes taken for parts of a glacial tongue	+	Usage of additional images for the period when the ice on the lake has already molten; usage of digital elevation models, in which this area will look such as an absolutely flat territory, contrasting with the sloped areas of the glacial tongue
Icings	+	Icing ice is usually lighter than glacial ice due to its cleanliness; such areas appear flat on digital elevation models

Figure 3. (A) Debris-covered part of Bolshoi Maashei glacier and the dead ice. 1—crevasses, marking the debris-covered part of the glacier; 2—thermokarst forms; 3—tunnels, both marking the dead ice areas: (B) debris-covered edge of glacier No. 138 (Aktru river basin); 4—the “smooth” debris surface; (C) seasonal snow cover and a large snowpatch (Maashei river valley); (D) the same area at the end of the ablation season.

Figure 3. (A) Debris-covered part of Bolshoi Maashei glacier and the dead ice. 1—crevasses, marking the debris-covered part of the glacier; 2—thermokarst forms; 3—tunnels, both marking the dead ice areas: (B) debris-covered edge of glacier No. 138 (Aktru river basin); 4—the “smooth” debris surface; (C) seasonal snow cover and a large snowpatch (Maashei river valley); (D) the same area at the end of the ablation season.

The delineation of recent glaciers and of the geomorphological forms that mark their positions in the maximum of the LIA has been completed manually. The minimum size of the glaciers to be mapped was 0.01 km².

The systematic error was defined as ±1 pixel (1.65 m for Geoeye-1, 1.8 for WorldView-2 and Corona, 10 m for Sentinel 2, 15 m for Landsat 7, 30 m for Landsat 5, and 40 m for Landsat 3). Thus, the area determination was calculated by a simple formula:

$$A\_er=n·m$$

(1)

where n is the number of pixels defining the perimeter of the glacier area, and m is the spatial resolution of the sensor bands applied expressed as an area of the pixel.

The percentage error of area determinations, A_rer, is given by:

$$A\_rer \left(\%\right)=A\_er/A\_\left(gl\right) 100$$

(2)

where A_gl is the area of the glacier.

To assess the dynamics of Lake Maashei, an analysis was made of aerial photographs of Corona for 1962 and satellite images of Landsat for 1980, 1989, 1993, 2000, 2004, and 2011, and data obtained during field work in September 2022 were used to evaluate the dynamics. The error assessment was made by the same algorithm as for the glaciers.

The main mistakes in determining the glacier outlines were associated with several factors that have already been described in our article, dedicated to the LIA glaciers and their shrinkage in Altai [26]. The most reliable way to avoid these mistakes is to use in situ observations whenever possible, but in cases where only remote sensing data are available, some indicator features can be used (Table 2).

We used a 30 m SRTM 1 Arc-Second Global DEM [28] to characterize our glacier outlines with parameters including mean, minimum, and maximum elevation ranges, and mean slope and aspect. Since our article also is dedicated to the reconstruction of the state of glaciers for the year 2000 and the maximum of the LIA, the use of SRTM DEM (around 2000) made it possible to better match the altitudes of the time points. The use of more modern DEM for glaciers in 2021 would not allow obtaining comparable data for all three time points. Such a comparison was possible provided that the data obtained were homogeneous. In addition, our experience of field work in the Altai mountains showed that the SRTM data most closely corresponded to the heights obtained in field studies using GPS.

Those parameters were determined automatically based on the DEM in the Global Mapper v.18.0 software (digitizer tool). Field data were used to verify the data obtained from remote sensing sources.

Since the equilibrium-line altitude (ELA) cannot always be detected on satellite images, especially for small glaciers, for the consistency of the obtained data on ELA, we used the Kurowski method [42,43] for its identification on glaciers excluding glaciers of valley type. In this method, the firn line altitude or ELA is calculated as the average altitude of the glacier:

$${\overline{z}}_f={\displaystyle \sum _i\frac{f_i{z}_i}{F}}$$

(3)

where ${\overline{z}}_f$ is the firn line altitude or ELA, f_i are the areas of the different altitudinal zones of the glacier, z_i are the average altitudes of these zones, and F is the total area of the glacier.

The Kurowski method is based on the assumptions that ablation and accumulation on a glacier change linearly with altitude and that the glacier is stationary. The assumption of linearity of changes in ablation with changes in altitude introduces a systematic error associated with the concave nature of the real curve of ablation versus altitude, due to which the ELA lies below the weighted average glacier elevation. At the same time, errors due to the assumption of linearity of changes in ablation with altitude and due to the assumption of glacier stationarity have the same signal during the period of glacier advance and opposite during the period of its retreat, that is, in the second case, they compensate each other, increasing the accuracy of the results. This exactly corresponds to the situation in our study, when glaciers were considered to be in the phase of intense retreat.

Therefore, the Kurowski method had high accuracy in relation to the recent glaciers of Altai and was also used to verify the values of ELA, found by remote sensing.

According to recent studies by Braithwaite [43], who tested the method for the 103 glaciers of different morphologies and from different regions of the world, there is a high correlation between balanced-budget ELA and Kurowski mean altitude, with a small mean difference of −36 m between the two altitudes with standard deviation ±56 m. The balanced-budget ELA is significantly lower (at a 95% confidence level) than the Kurowski mean altitude for outlet and valley glaciers, and not significantly lower for mountain glaciers.

In our recent work, we tested the Kurowski method for 26 mountain glaciers of the Shapshalsky Center (cirque and hanging glaciers with areas less than 0.6 km²) [26,44]. Calculations using this method gave results close to the ELA position obtained by remote sensing (the average value of the difference was −8 m). Since Braithwaite [43] indicated that the accuracy of the method was noticeably reduced for large glaciers, we performed an additional study for the glaciers of the Tavan Bogd massif: 26 glaciers in the area range of 3.3–23.1 km² [26]. In that case, the calculated ELA values by the Kurowski method as of 4 points in time and the data obtained from the corresponding satellite images were considered. We obtained an average difference of 77 m, which means general overestimating of the ELA altitude by the Kurowski method occurred. Taking into account the vertical span of the glaciers studied, which reached about 1.5 km, such an error can be considered acceptable.

Nevertheless, we conducted an additional study on the applicability of the Kurovsky method for large glaciers using the example of the valley glaciers (35 glaciers) of the North-Chuya ridge (

Table 3. The difference between the ELA calculated by the Kurowski method (ELA_k) and obtained from satellite imagery (ELA₂₀₂₁) for the valley and the characteristics of the glaciers.

No (Name)	ELAk-ELA2021	Area, km	Length (L), km	Vertical Extent ΔZ (m)
3 (Jelo)	25	6.34	3.72	870
15	18	0.36	1.49	597
19 (Malyi Aktru)	171	2.54	3.63	1223
22 (Pravyi Aktru)	187	4.25	4.18	1191
25 (Levyi Aktru)	205	4.93	5.34	1271
30 Yan-Karasu	101	1.09	2.96	1155
31 (Malyi Korumdu)	0	0.36	2.15	1036
32 (Korumdu)	180	5.52	4.79	1607
34 (Kurkurek)	257	1.96	3.36	1407
39	261	1.97	2.84	1115
42 (Pravyi Maashei)	218	2.89	3.57	1195
46	106	4.84	4.43	1157
47 (Bolshoi Maashei)	390	6.22	7.82	1882
48	363	2.12	2.85	1476
51 (Levyi Maashei)	283	3.46	2.87	1353
69	116	0.69	1.56	742
71	76	1.65	2.77	1001
102	34	1.81	2.84	868
104	139	2.44	2.50	1083
109	−29	0.67	1.60	446
111	248	2.13	3.42	1202
116	−43	1.04	2.84	897
143	4	0.66	1.84	650
144	−6	0.33	1.24	638
181	44	0.54	1.56	395
191	−51	0.62	1.65	599
199	−40	0.85	1.45	454
200	13	1.14	2.20	653
202	−32	0.67	1.57	359
203	102	2.18	3.12	977
210	219	4.22	4.95	1213
212	292	2.37	3.06	1370
213	233	0,82	2.17	1166
214	176	3,11	3.12	957
215	221	2,64	3.39	965

).

The average difference between the ELA calculated by the Kurowski method and the ELA obtained from the satellite image (128 m) was generally higher than the values that we had obtained earlier or the Tavan Bogd massif. At the same time, a correlation was found between the ELA_k-ELA₂₀₂₁ values and the glacier area (0.53) and the glacier length (0.63), and the highest correlation was found with the vertical extent of the glacier (difference between the altitude of the highest and the lowest points of the glaciers) −0.73.

Based on the results of the comparison (Table 3), for valley glaciers, we included in the glacier inventories of 2021 and 2000 the results obtained directly from the interpretation of images; for glaciers of other morphological types, the position of the ELA was obtained by the Kurowsky method. Unfortunately, for the glaciers reconstructed for the LIA maximum, it was not possible to use satellite imagery. Therefore, for the valley glaciers of the LIA we used the correction calculated on the basis of the linear dependence between the vertical diapasone of the glaciers ΔZ(m) and ΔELA = ELA_k − ELA_2021 (Table 3).

$$\Delta ELA=0.2945\ast \Delta Z-168$$

(4)

3.3. Paleo Reconstructions

The reconstruction of the LIA extent of the glaciers was completed on the basis of the geomorphological methods. The LIA moraines were mapped using satellite imagery, GPS tracking of the lateral and terminal moraines, and visual in situ observations. We used the method of ground-based route interpretation that included descriptions, measurements, and photography in the reference areas. For object recognition, we used the visual interpretation method according to the reference standards [45,46]. In our previous works [26,47,48], we compiled the interpretation standards on the basis of ground-based observations, indicating the following characteristics: characteristic images of objects on the terrain, in the aerial photograph, in the satellite image, in the topographic map in the DEM; distinctive characteristics of objects; and methods of transferring objects to the map. Working with the satellite imagery and DEM, we used the criteria similar to those suggested in [49]: “Identification criteria include shadowing due to changes in topography (relative relief) and changes in color due to changes in soil, soil moisture, and vegetation cover. Associated landforms, such as deflected abandoned meltwater channels, are also useful in delineating the break-of-slope of these features”.

According to our interpretation standards, the diagnostic features of LIA moraines are as follows:

-

Steep fronts and relatively large thickness;

-

Glacial ice cores that are sometimes exposed by thermokarst processes and form thermokarst depressions, landslides, and thermoerosional forms, such as the tunnels in the dead ice, that have been described in Table 2 and shown in Figure 4;

-

Position adjacent to modern glaciers. This criterion cannot always be used because some of the glaciers of the LIA completely disappeared. However, in some cases, we can find some of the glaciers that have already disappeared in the Landsat images of 2000 or the Corona images of 1968 and 1962;

-

The low degree of vegetation cover of the LIA moraines, which in multispectral images was expressed by a grey or brown color of moraines in sharp contrast with the greenish color of the surrounding subalpine meadows and tundra and the moraines of earlier glacier advances, so they are easily identified. It should be noted that this diagnostic feature for the North-Chuya ridge needs to be clarified. In the relatively humid conditions of the North-Chuya ridge, the overgrowth rate of moraines is higher; therefore, the lower parts of moraine complexes can be grass-covered or even covered with trees. At the same time, for the lateral moraines of the Little Ice Age, this overgrowing has not yet manifested itself; therefore, in such cases, it is necessary to track the moraine lines towards the glacier front.

Sometimes the similarity of the LIA moraines to rock glaciers leads to an overestimation of the glacier area. The ways to solve this problem have been described in our previous article [26]: in the case of talus rock glaciers [50], they are located below talus slopes devoid of glacial exaration forms such as cirques or corries, and debris rock glaciers develop below the LIA moraines or overlap them and move further down the valley. In cases of overlapping of the LIA moraines by the debris rock glaciers, the lateral moraines are still overlooking them, and it is usually possible to extrapolate them as they approach each other down the valley to find the lowest point of the LIA glacier. When the lateral moraines are not clearly visible, a combination of a DEM and a satellite image (3D model) helps to solve the problem.

Besides moraine mapping when reconstructing the glaciers, we used the clear boundaries between unweathered freshly glacially eroded areas and weathered, vegetated areas, in particular fresh glacial erosion marks on the walls of the cirques and troughs in the upper parts of the valleys, corresponding to the LIA moraines located further below the valley.

The diagnostics of the hanging glaciers were carried out on nival niches with sharp outlines, indicating the recent degradation of the glaciers that formed them. As a rule, such niches are marked with modern snow patches. For some of these glaciers, the degradation process has occurred recently and was detected by satellite images and aerial photographs (Figure 5), making it possible to use the nival niches that remained in their place as interpretation standards.

As we have shown above, the observations of the glaciers of Altai in the 19th century led to the conclusion that the climax of the glaciers took place between 1835 and 1880, most likely around 1850 [24,35]. However, it is not obvious whether this advance was maximal during the LIA.

The timing of the maximal LIA glacier advance in Altai is interpreted differently by different researchers: in the end of the 15th century [51], in the 17th century [52], and in the beginning [53] or in the middle of the 19th century [27,54]. Most researchers compare the LIA maximum in Altai to the Aktru stage in the 17–19th centuries [24,55]. According to Okishev [24], during that stage there was a two-fold glacial advance of equal scale in the early 17th and mid-19th centuries. This is in good agreement with our earlier estimations of 1850 as the starting point for the general retreat of the Altai glaciers [26], based on the mass balance index calculations and the estimations of the time of reaction of the termini positions of the valley glaciers of Mongun-Taiga and Tavan Bogd.

3.4. Estimation of the Thickness and Volume of Glaciers

Methods for estimating the thickness and volume of glaciers in the framework of our study can be divided into two groups: (a) empirical area-volume scaling obtained as a result of radar sounding; and (b) glaciological modeling of ice thickness, based on the physical laws of their movement, data on boundaries, and surface topography.

a.

Empirical area-volume scaling.

Based on the assumption that the transverse profile of the valley glacier has a parabolic shape, V.N. Erasov [56] suggested that its volume should be proportional to the area in the form of expression (5):

$$V=k{S}^p$$

(5)

where S is the area, k = 1.5, and the coefficient p is related to the morphometric and morphological characteristics of a particular glacier.

S.A. Nikitin [33] obtained regional dependences between the area and the volume of glaciers on three ridges of the Russian Altai: Katunskiy, North-Chuya, and South-Chuya. Based on the same data set, Macheret [57] derived an empirical equation for flat-topped glaciers.

Table 4. Dependence between the area (S) and volume (V) of Altai glaciers as a function V = kS^p.

Morphological Type	K	p	Source
Valley (n = 46; r = 0.94)	0.0435	1.165	[33]
Cirque-valley (n = 36; r = 0.89)	0.0464	1.088
Cirque and cirque-hanging (n = 29; r = 0.91)	0.0487	1.244
Flat summit (n = 5; r = 0.989)	0.044	0.89	[57]

presents the regional Altai values of the coefficients for various morphological types of glaciers that were used in our study.

Using expression (5), we also determined the average thickness of the glacier:

$$H\_av=V/S$$

(6)

where V is the volume of the glacier and S is its area.

The maximum thickness of the glaciers was also determined using the formula proposed in [58]:

$$H\_max=77.5S^{0.245}$$

(7)

where S is the area of the glacier.

b.

Ice thickness modeling.

Modeling the thickness of glaciers, based on the physical laws of their movement, data on their boundaries, and their surface topography, presented in the form of digital elevation models (DEM), is currently the most promising direction for assessing the volume of glaciers’ potential water reserves [59,60].

The simple and easy-to-use model glacier bed topography (hereafter referred to as GlabTop) [61,62] and its modifications (the result of development by the authors of the original model-GlabTop2 [62]) have been successfully used in many regions of the world. The model assumes that the shear stress on the bed (τ) along the central line of the glacier is constant [63], the flow of the glacier is laminar, and the contribution of modern glaciers to the formation of their bed is small compared to the existing subglacial relief, and that the underlying surface is smooth [64]. Thus, the thickness of ice is largely determined by the slope of its surface (the steeper the surface, the thinner the ice, and vice versa). Ice thickness along the glacier center lines is estimated from the equation:

$$h=\tau /\left(fg\rho sin\alpha \right)$$

(8)

where h is the ice thickness; τ is the shear stress on the bed; f is the cross-sectional shape factor of the glacier; ρ is the ice density; g is the free fall acceleration; α is the angle of inclination of the surface along the central line of the glacier.

The main difference between GlabTop2 and GlabTop is that the slope of the surface is calculated not along the glacier center lines, but as the average slope of the surface of all pixels in the raster within a certain buffer. This allows one to fully automate the calculations [65].

GlabTop2-py [github.com; pypi.org] is a Python package that computes the ice thickness distribution for all glaciers in a simulated area. GlabTop2-py uses the Python 3.8 [python.org] and PCRaster [pcraster.geo.uu.nl] features. The model is fully based on the concepts described in [65].

The simulation result is a raster image of the spatial distribution of ice thickness. With the help of standard statistical GIS tools, the average and maximum thicknesses of the glacier were determined.

3.5. Glacioclimatic Calculations

To determine the climatic changes that determined the nature of the shrinkage of the Bolshoi Maashei glacier, we calculated the mass balance index according to the method of G.E. Glazyrin [66], in which the calculation requires data on air temperature and precipitation at the base meteorological stations, as well as the ELA:

$$I\_b=c-a$$

(9)

where I_b is the mass balance index; c is the accumulation; ɑ is the ablation; all calculations of these quantities were performed in mm w.e.

To calculate ablation at any point of the glacier, it is necessary to calculate the average summer temperature t_i at the corresponding altitude z_i using the vertical temperature gradient G_t, the altitude of the weather station z_m, and the average summer temperature at the meteorologic station.

For calculations, we used data from the Kara-Turek weather station, located approximately 82 km west of the glacier at an altitude of 2600 m a.s.l. The vertical temperature gradient was calculated using a pair of Akkem-Kara-Turek meteorological stations (the distance between these stations is approximately 14 km).

The temperature jump between the glacial and non-glacial surface also should be taken into account. According to [67], for glaciers in the range of areas from 4 to 10 km², the temperature jump is −1.5 °C.

As a result, we obtained:

$$t\_i=t\_m-G\_t \left(z\_i-z\_m \right)-\Delta t$$

(10)

After calculating t_i, the ablation was determined. For this, we used the ablation formula according to the improved variant of formula by A.N. Krenke and V.G. Khodakov [68], supplemented by the regional exposure coefficient k [69]:

$$\mathsf{\alpha }=k1.33\left(t\left(z\_\mathrm{i} \right)+9.66\right)^2.85$$

(11)

For the glaciers with the northern aspect, this coefficient was 0.82.

Accumulation was calculated by the formula

$$\mathrm{c}=KP$$

(12)

where K is the concentration coefficient; P is the annual rainfall.

To determine the annual precipitation by extrapolation from the Kara-Turek weather station, it is necessary to know the vertical pluviometric gradient, which was obtained from a pair of Kara-Turek-Akkem weather stations: its average long-term value was 11.5 mm/100 m. The concentration coefficient was calculated using the glacier dependence from the area for Altai obtained by V.P. Galakhov [69]:

$$K=1.3438\times 9.17S^\left(-0.141\right)$$

(13)

where S is the area of the glacier.

For the Maashei glacier, the K value turned out to be close to 1.

4. Results

4.1. Reduction of the Glaciers of the North-Chuya Ridge after the Maximum of the LIA

To assess the reduction of glaciers from the maximum of the LIA, glacier inventories were compiled for the maximum of the LIA in 2000 and 2021.

The total systematic error for the 2021 glacier inventory was 8.1%, and for 2000 it was 10.0%. To evaluate the subjective ‘cartographer’ error, we took a sample of 20 glaciers. Within the sample, the distribution of glaciers by area corresponded to the distribution by area of glaciers within the entire data set. Sample glaciers were remapped. The error was determined by a comparison of the remapping results with the glacier areas obtained over the course of cataloguing. The average error for the entire sample was 4.9%.

For the maximum of the LIA, we reconstructed a total of 400 glaciers with a total area of 304.90 ± 23.49 km², and the area-weighted mean ELA was 2909 m a.s.l. (

Table 5. LIA glaciers of the North-Chuya ridge. Designations in the table: N—glacier number according to this inventory; S—total areas of the glaciers of the river basin, km²; Z₁—the altitude of the lowest point of the glaciers, m a.s.l.; Z₂—the altitude of the upper point of the glacier, m a.s.l.; E—average aspect of the glacier.

River Basin	N	S	Z1	Z2	ELA	E
Jelo	12	17.73 ± 0.73	2512	3695	3050	SE
Tete	15	8.92 ± 0.53	2641	3551	3099	NNE
Aktru	20	45.02 ± 2.25	2189	3970	2840	NE
Maashei	19	54.83 ± 2.56	2196	4110	2902	N
Shavly	43	38.23 ± 2.56	2157	3800	2823	NW
Yungur	138	46,67 ± 3,98	2219	3525	2902	WNW
Argut	10	5.09 ± 0.39	2259	3342	2884	SW
Karagem	143	88.41 ± 5.70	2262	4113	2945	SSE

; Table S1, Supplementary Materials). The total volume of the glaciers was approximately 15.5 km³ according to the results of calculations using the formulas of Nikitin and Macheret. The largest of the glaciers reached an area of 17.49 km² (Bolshoi Maashei), the lowest point of the glacial expansion was 2157 m a.s.l. (Shavly river basin), and the lowest ELA was 2341 m a.s.l. (Yungur river basin). The Karagem river basin was the most glaciated: 143 glaciers with a total area of 88.41 km², although the glaciers of the Maashei and Aktru basins were two–three times larger on average, though less numerous.

By 2000, the total area of the glaciers decreased to 140.24 km² (54% decrease, 1,1 km² (0.36%) per year), and the number of the glaciers decreased to 224 (44% decrease) (

Table 6. Glaciers of the North-Chuya ridge in 2000. Abbreviations are given in Table 5.

River Basin	N	S	Z1	Z2	ELA	E
Jelo	11	10.78 ± 0.97	2694	3688	3150	SE
Tete	7	2.72 ± 0.41	2804	3551	3158	NNE
Aktru	22	28.70 ± 3.20	2255	3970	3093	NE
Maashei	37	34.32 ± 3.32	2205	4110	3037	N
Shavly	42	20.51 ± 2.74	2266	3803	3009	NW
Yungur		10.75 ± 1.76	2581	3525	3050	WNW
Karagem	54	32.45 ± 3.80	2432	4113	3111	SSE

; Table S2, Supplementary Materials). In total, 273 glaciers disappeared completely. The volume decreased to approximately 6.9 km³ (55.5% decrease) according to the calculations by the Nikitin/Macheret formulas, and to approximately 6.6 km³ by the GlabTop method. The area-weighted mean ELA in 2000 reached 3074 m a.s.l. (165 m rise since the LIA maximum). The lowest point of the glaciers reached 2205 m (an increase of 48 m), and the lowest ELA was 2670 m a.s.l. (269 m rise). The Karagem river basin still had the largest total area of the glaciers, though it became almost equal to that of the Maashei river basin (only 0.1 km² difference).

In the period of 2000–2021, the total area of the glaciers decreased to 120.02 km² (

Table 7. Glaciers of North-Chuya ridge in 2021. Abbreviations are given in Table 5.

River Basin	N	S	Z1	Z2	ELA	E
Jelo	10	8.47 ± 0.52	2818	3688	3179	SE
Tete	6	1.62 ± 0.13	2880	3551	3208	SE
Aktru	21	25.34 ± 1.99	2357	3970	3134	NE
Maashei	36	30.79 ± 2.06	2225	4110	3081	N
Shavly	44	17.31 ± 1.65	2440	3803	3059	NW
Yungur	55	8.17 ± 1,07	2641	3525	3075	WNW
Karagem	52	28.31 ± 2.29	2525	4113	3160	SSE

; Table S3, Supplementary Materials). The loss of area in 2000–2021 made 7% of the total area in the LIA maximum, 0.96 km² decrease per year on average, and 10 glaciers disappeared. The total volume of the glaciers decreased to 5.8 km³ (16% decrease since 2000 and 65% decrease since the LIA maximum) by the Nikitin/Macheret formulas, or to 5.6 km³ (15% decrease since 2000) by the GlabTop method. The marked slowdown in the reduction of the glacier area was largely due to the degradation of a large number of small glaciers at the stage of 1850–1920; by 2000, with a decrease in the number of small glaciers, the potential for a rapid reduction in the area of glaciers decreased. Despite the decrease in the absolute values of the reduction in the area of glaciers, their relative reduction accelerated: 0.36% decrease of the area per year in 1850–2000; 0.69% decrease of the area per year in 2000–2021.

The area-weighted mean ELA in 2021 was 3116 m a.s.l. (42 m higher than in 2000 and 207 m higher than in the LIA). The lowest point of the glaciers reached 2225 m a.s.l. (a rise of 20 m since 2000), and the lowest ELA was 2718 m a.s.l. (48 m rise since 2000). The Maashei river basin had the largest total area of glaciers (more than 25% of the total for the North-Chuya ridge).

An analysis of the reduction of glaciers on the North-Chuya range from the maximum of the LIA to 2021 (

Table 8. Reduction of glaciers of the North-Chuya range from the maximum of the LIA to 2021 by river basins. ΔS km (%): changes of the total glacier area, km (% to the LIA area); ΔN: changes of the number of glaciers (% to the LIA); N_d: the number of the glaciers that disappeared (% to the LIA).

River Basin	S	E	ΔS km (%)		ΔN		Nd
			1850–2000	2000–2021	1850–2000	2000–2021	1850–2000	2000–2021
Jelo	1.48	SE	−6.95 (−39%)	−2.31 (−13%)	−1 (8%)	−1 (8%)	4 (33%)	1 (8%)
Tete	0.64	NNE	−6.2 (−70%)	−1.1 (−12%)	−8 (53%)	−1 (7%)	10 (67%)	1 (7)
Aktru	2.25	NE	−16.32 (−36%)	−3.36 (−7%)	+2 (10%)	−1 (5%)	8 (40%)	1 (5%)
Maashei	2.89	N	−20.51 (−37%)	−3.53 (−6%)	+18 (95%)	−1 (5%)	5 (26%)	2 (11%)
Shavly	0.89	NW	−17.72 (−46%)	−3.20 (−8%)	−1 (2%)	−6 (14%)	19 (44%)	2 (5%)
Yungur	0.34	WNW	−35.92 (−77%)	−2.58 (−6%)	−87 (63%)	+4 (3%)	108 (78%)	1 (<1%)
Argut	0.51	SW	−5.09 (−100%)	-	−10 (100%)	-	10 (100%)	-
Karagem	0.62	SSE	−55.96 (−63%)	−4.14 (−5%)	−89 (62%)	−2 (1%)	109 (76%)	2 (1%)

, Figure 6, Figure 7, Figure 8 and Figure 9) by river basins revealed that until 2000, small glaciers were predominantly reduced and degraded, respectively, and the area of glaciers in river basins with the smallest average area of the glaciers (Argut, Yungur, Tete, and Karagem) underwent the greatest decrease, mostly due to the disappearance of a large number of glaciers. However, in the period 2000–2021, such a pattern did not manifest itself, and the reduction across the basins occurred relatively evenly, which probably reflects an increase in the rate of reduction of relatively large valley glaciers that previously showed resistance to climate change due to their greater inertia.

Glaciers of medium and small sizes have undergone severe degradation, although this is most obvious for medium-sized glaciers (the cirque-valley type). While degrading, the cirque-valley glaciers were transformed into cirque glaciers, which compensated for the degradation of this type of glacier. At the same time, the degradation of valley glaciers and their transition to cirque-valley glaciers were manifested to a lesser extent. Accordingly, it was the cirque-valley glaciers that underwent the greatest area reduction (Figure 10). The relative stability of valley glaciers manifested itself in an increase in their share in the overall structure of glaciation, against the background of a relatively greater reduction in small- and medium-sized glaciers.

The aspect distribution of the glaciers of the North-Chuya range has changed little (Figure 11). The glaciers of the north-western aspect, which had the smallest area at the maximum of the LIA, underwent the least reduction. The predominance of glaciers with the northern aspect became more pronounced because the glaciers of the north slopes decreased to a lesser extent compared to the glaciers of other slopes. A possible reason for the different glacier reduction of the shaded northern slopes and the more insolated slopes of other aspects could be the growth of insolation difference between them. This could happen as a result of the general increase in insolation against the background of a decrease in precipitation. Such an aridization trend has been traced at least since the 1960s, when systematic observations began (Figure 12).

Unfortunately, weather stations on the territory of the Altai mountainous country cover only the period from the 1960s, which does not allow us to assess the change in climatic conditions from the maximum of the LIA. The nearest meteorological station with a series of observations partially covering the LIA is Barnaul. The analysis showed a high correlation (correlation coefficient of more than 0.7) of the values of the average summer temperatures of the Barnaul and Kara-Turek meteorological stations for the period 1961–2021, respectively. Consequently, we used linear dependence obtained between these data:

$$T_K=0.82\ast T_B-8.94$$

(14)

where T_K_ is the average summer temperature for the Kara-Turek station, T_B.

Using (14), we extended the data series on the average summer temperatures for the Kara-Turek station until 1838 (Figure 13), which made it possible to characterize changes in thermal conditions in the study area in the most general way.

According to the data obtained, the trend toward an increase in summer temperatures manifested itself approximately until the early 1920s, and it was most pronounced during the last 30–40 years of this time interval. From the beginning of the 1920s to approximately 1970, the trend in summer temperature changes was negative, which was a probable cause of the slowdown in the reduction of some large glaciers (for example, the Bolshoi Maashei glacier; see the next section).

An analysis of data from the nearest meteorological stations Kara-Turek (Figure 10), Kosh-Agach, and Akkem shows that against the general background of warming resumed since the 1970s, a sharp increase (by 1.5–3 °C) in average summer temperatures in the 1990s was visible in all records, and from the early 2000s there has been some stabilization at a high level. This ‘thermal shock’ undoubtedly caused an acceleration in the degradation of glaciers, primarily small ones. For larger valley glaciers, this event has a delayed effect, which began to fully manifest itself only in the last decade, which we will consider in more detail below using the Maashei glacier as an example.

The fact is that no acceleration in the reduction of the total area of glaciers in the last 20 years has generally been noted, probably due to a certain inertia of large glaciers, while their relative share (as can be seen in the example of valley glaciers) has increased greatly due to the degradation and disappearance of small glaciers in the previous time interval.

4.2. Dynamics of the Bolshoi Maashei Glacier

At the time of the maximum of the Little Ice Age, according to our reconstruction, the Bolshoi Maashei glacier had an area of 17.49 km².

The glacier consisted of three streams (Figure 14), merging at altitudes of about 2650–2700 m a.s.l. The largest in the vertical range was the central stream, originating in the cirque of eastern exposure right under the highest point of the North-Chuya ridge (4177.7 m a.s.l. Mount Maashei-Bash). The length of the glacier along this stream was 8.9 km. The width of the stream was 850–900 m, narrowing at altitudes from 3000 to 2900 m a.s.l. to about 300 m. Further downslope, it occupied another cirque with a bottom level of approximately 2700–2750 m a.s.l. The glacier along the eastern stream had approximately the same length, while it had a large width (from 850 to 1500 m); however, at the same time, it began relatively low (at an altitude of about 3750 m a.s.l.). The predominant direction of this stream, before its confluence with the central one, was west-north-west. The shortest western stream began at a level of about 4090 m a.s.l. and had a length of about 2.7 km to the point of confluence with the central stream. The ice here, descending the steep slope of the north-eastern exposure, completely occupied the cirque of the north-eastern exposure with bottom marks of about 2900 m a.s.l. and, after leaving it, merged with the central stream.

After the confluence, all flows formed a single glacial tongue that spanned about 3.9 km. The lowest point of the glacier was approximately at an altitude of 2160 m a.s.l.

After the start of observations, Tronov and Okishev noted low rates of glacier shrinkage. At the same time, Tronov noted the formation of a small moraine rampart in the late 1920s [35], and P.A. Okishev noticed another shaft around 1944–1947. The remains of these swells were observed and are currently mostly closer to the eastern side of the valley, and they are completely eroded along the front. On the western side, we noted the remains of a moraine rampart dating back to 1947, no more than 1 m high (Figure 15); however, a similar state of it was noted here by P. Okishev in 1962 [32].

According to the USSR Glacier Catalog for the 1960s [29], the glacier remained intact, and its area was 16.0 km². At the same time, our analysis of the Corona images of 1962 and 1968 showed that the eastern flow of the glacier had already lost contact with the main glacier by 1962 (Figure 16). M.V. Tronov, who observed the glacier in 1937–1939 noted that the glacier was already retreating and this stream, forming a wide and flat tongue, was largely isolated, and only its left part, which had additional snow avalanche accumulation, merged with the middle stream [70] (p. 147).

According to our data, in 1968, the area of the Bolshoi Maashei glacier was 9.17 km², and the eastern stream that separated from it became an independent glacier (hereinafter referred to as Eastern Maashei) with an area of approximately 6.09 km². Thus, the reduction of the entire glacial system from the LIA maximum amounted to 2.23 km², and approximately half of this value was not associated with the reduction of the glacier tongue, but with the retreat and separation of its eastern flow.

The next important event in the dynamics of the glacier occurred between 1993 and 2000, when the glacier finally lost contact with its western flow. In 1993, according to our data, the area of the main glacier was 8.76 km² (average reduction 0.18%/year for 1968–1993), the Eastern Maashei glacier had an area of 5.93 km² (reduction 0.10% per year), and the total reduction of the glacier system was 0.57 km² (0.15%/year on average).

By 2000, the area of the detached western flow, which became an independent glacier (hereinafter West Maashei) was 2.15 km², the Bolshoi Maashei glacier decreased to 6.45 km² (the area of reduction for the initially unified glacier here was 0.16 km², 0.26%/year on average), and the Eastern Maashei glacier was 5.21 km² (1.7% per year on average). There was a sharp acceleration in the loss of area by the entire glacial system (0.86% per year), and it is noteworthy that this mainly happened not due to the degradation of glacier tongues, but due to the collapse of the Bolshoi Maashei glacier and the degradation of parts of the Eastern Maashei glacier in the altitudinal interval of 2700–3400 m a.s.l., where accumulation had been carried out mainly due to avalanches. This fact is probably associated with a sharp reduction in slope snow cover of the slopes; a similar reduction in glaciers above glacial tongues for the same years was observed earlier for the southeastern Altai [25].

By 2022, the area reached 2.12 km² for the Western Maashei glacier (0.06% reduction per year in 2000–2021), 6.25 km² for the Bolshoi Maashei glacier (0.14% reduction per year), and 4.84 km² for the Eastern Maashei glacier (0.32% reduction per year). Therefore, the area of the entire glacial system has decreased by 0.2% per year for the last 22 years. Apparently, such a slowdown in relation to the previous time interval is due to the fact that in the previous stages, the glacial system disintegrated by detaching glacial flows, which gave a sharp reduction in the area of glaciers; now, there is only a retreat of glacial tongues, which gives a slower area decrease of the glacial system. In addition, the Western Maashei glacier has retreated inside the cirque, where the steepness of the surrounding slopes provides shading and the additional accumulation of avalanche snow on the tongue, which has slowed its decline. The relatively high rate of reduction of the Eastern Maashei glacier is apparently associated with its lowest average hypsometric position and the greater width of the enclosing trough, which contribute to the accelerated degradation of the glacier at the present stage.

As of 2022, the tongue of the Bolshoi Maashei glacier extends northward for about 3.9 km. In the upper part, its width reaches 900 m, gradually narrowing to 500–600 m. Both the eastern and western edges of the glacial tongue are debris-covered. The two marginal bands of debris-covered ice reach in some parts up to 300 m wide. They have not yet lost their movement, but they are already clearly showing a tendency to turn into dead ice and to get separated from the more dynamic central part of the glacial flow by water streams. Obviously, the moraine cover significantly reduces the ablation of contaminated areas of the glacier; due to this, they rise 10–15 m above the open ice strip in the lower part of the glacier, which enhances the feeling of a glacial stream flowing in the stone-ice shores.

At the end of the 2022 ablation season, the edge of the Bolshoi Maashei glacier was located at an absolute altitude of 2225 m a.s.l. (Figure 17). The central part of the glacier, slightly contaminated with moraine material, descends the lowest. The Maashei river breaks out directly from under the edge of the glacier. The front of the glacier has the character of a steep forehead; visually, its thickness is up to 15 m. The western part of the tongue ends at approximately the same elevation.

The ELA for 2021, determined from the Sentinel 202109/08 satellite image, averaged 2890 m.

When analyzing data on the reduction of glaciers, we had to take into account a rather large measurement error associated with the low resolution of images from the 1980s–1990s. Therefore, in our opinion, it is more significant to consider longer time intervals for which the relative error in determining the retreat of the glacier is small (

Table 9. Retreat of the tongue of the Bolshoi Maashei glacier.

Period	Retreat for a Given Time Period	Average Retreat Rate, m/Year	Averaging Period	Retreat for a Given Time Period	Average Retreat Rate, m/Year
1850–1924	320 ± 5 1,2	4.3	1850–1924	320 ± 5 1,2	4.3
1924–1932	49 2	6.1	1924–1962	244 2,3	6.4
1932–1937	35 2	7.0
1937–1952	75 2,3	5.0
1952–1962	85 3	8.5
1962–1968	24 ± 1 3,5	4.0	1962–1989	87 ± 30 5	3,2 ± 1
1968–1975	28 ± 1 5,3	4.0
1975–1980	11 3,4	2.2
1980–1989	24 ± 33 4,5	2.7
1989–1993	35 ± 60 5	8.8	1989–2010	148 ± 45 5	7.1 ± 2.0
1993–2000	40 ± 45 5	5.7
2000–2004	16 ± 30 5	4.0
2004–2007	30 ± 30 5	10.0
2007–2010	27 ± 30 5	9.0
2010–2013	59 ± 17 5	19.8	2010–2022	208 ± 18 5	14.0 ± 1.5
2013–2016	42 ± 4 5	14.0
2016–2019	21 ± 12 5	7.0
2019–2022	86 ± 13 5	28.6

Notes. ¹—according to geomorphological data; ²—field data, M.V. Tronov [70]; ³—field data, P.A. Okishev [24,32]; ⁴—field data, R.M. Mukhametov; ⁵—based on the interpretation of satellite images.

). The duration of the first interval (1850–1924) was determined approximately, since the exact date of the culmination of the Bolshoi Maashei glacier in the LIA is unknown, although based on its very large size for Altai, one can assume some inertia of its behavior, and the beginning of the retreat can be attributed only to 1850. The retreat values obtained by us for 1850–1924 are in good agreement with the low rates of retreat of the glacier established for it for all subsequent periods, with the exception of the last 12 years (Table 9), when it accelerated retreat under the influence of extreme warming. The next interval (1924–1962) can be characterized as a time of slightly higher glacial retreat rates on average, but these rates were uneven, as evidenced by the formation of moraines in the late 1920s and in 1944–1947. Probably, during these years, the rate of glacier retreat approached 0.

The period of decrease in average glacier retreat rates in the period 1962–1989 was probably associated with a decrease in average summer temperatures that took place from the early 1920s to approximately the end of the 1960s (Figure 12).

The subsequent general trend toward a gradual acceleration in the retreat of the glacier edge can be associated both with the consequences of the collapse of the previously unified glacier and a decrease in the flow of ice to its tongue due to the separation of lateral glacial flows and with the progressive deterioration of the climatic conditions for the existence of glaciers. It should be noted that it is difficult to tie the first process to specific time points, since the reduction in ice inflow from lateral flows probably occurred gradually, and the effect of this process was extended over time.

The results of calculations of the mass balance index changes according to (9)–(13) (Figure 18) show that from the beginning of the 1960s to the mid-1980s, there was a weakly positive trend in the change in the mass balance; from the mid-1980s, the trend changed to a pronounced negative one, which continues to this day.

Interestingly, this negative trend was associated with both an increase in summer temperature and a decrease in annual precipitation (Figure 9); moreover, the temperature increase was very pronounced in the 1990s and slowed down in the last 20 years (a similar pattern is also typical for weather stations in other regions of Altai), and the amount of precipitation continues to decrease.

The results of the calculations of the mass balance index and data on temperature changes suggest that the first impulse to accelerate the retreat of the glacier since the late 1980s was associated with an increase in summer temperatures and a sharp increase in the melting of the tongue part of the glacier (according to our calculations, melting in the lower part of the tongue increased from the end of the 1980s to the end of the 1990s by 1 m w.e., and in 1998, the calculated melting reached 5.1 m w. e.). It is likely that the second impulse to accelerate the retreat of the glacier around 2010 was when less-thick ice, which formed in the accumulation zone during the years of mass balance deficit in the late 1980s and 1990s, approached the glacial front.

The accelerated melting and retreat of the Bolshoi Maashei glacier has had a significant impact on the surrounding landscapes. The Bolshoi Maashei glacier is one of the largest and longest glaciers in the Russian Altai with almost 4 km located below the ELA; therefore, its melting is one of the significant factors in the occurrence of mudflows in the Maashei river valley.

The combination of large ablation values and high amounts of liquid precipitation is especially favorable for high runoff. In addition, ablation can increase due to heavy rainfalls, which adds up to the already high influx of additional water into the valley. In the presence of dammed lakes in the hydrological system of the valleys, this flow is regulated; however, in the event of the destruction of the dam and/or the overflow of the lake basin, especially destructive mudflows occur. This is facilitated by the presence of loose material easily transported by water on the bottoms of glacial valleys.

In the Maashei valley, there was a relatively large lake fed by glacier runoff at a distance of about 5 km from the front of the Maashei glacier (Figure 19). It was dammed by a rock glacier, sliding down from the western slope, and blocked almost the entire valley along with a talus that joined with it at the foot of the eastern slope of the valley.

According to our interpretation of the satellite images of different years between 1993 and 2012 and our in situ geodetic survey, based on the geomorphic evidence of high-level stands (Figure 20), the parameters of the lake varied during the seasons and in different years from complete drainage to the maximum before its outburst in 2012, when its length was 1480 m, width 425 m, the maximal area 259,360 m², and water volume was about 1,212,210 m³. The maximum depth of the lake according to the survey results was 7.5 m, with an average of 4.7 M. This is more than the average depth that had been reported earlier: 3–3.5 m [71].

In the summer period of 2012, the overflow of the lake basin, the outburst of the lake, and the passage of a catastrophic debris flow along the Maashei valley were caused by the coincidence of high values of two meteorological factors at once: temperature and precipitation. The average monthly temperatures in June and July that year exceeded the long-term average by 3.4° and 1.2°, respectively (Kara-Turek station). The ablation period of 2012 was one of the most negative for the glacier mass balance. According to our calculations, the melting for this period was estimated to be about 2000 mm w.e. at the ELA and 4800 mm w.e. at the glacial terminus. The amount of precipitation according to the Kosh-Agach weather station, located in the driest region of the Russian Altai, for the second decade of July amounted to 321% of the average multiyear norm, while only from July 12 to 15, 41.7 mm fell (link to the site of the Kosh-Agach weather station), with an average long-term value of about 120–150 mm/year (link to the site of the weather station or the work of Rusanov 1961). In the Kara-Turek station, the daily precipitation in the period before the outburst was the following: July 12—12.6 mm, July 13—13.6 mm, July 14—36.7 mm [72].

The superimposition of peak values of atmospheric precipitation and glacier ablation led to the active removal of detrital material from glacial valleys. This process was especially significant in the Katunsky, South-Chuya, and North-Chuya ranges. On 15 July in the Maashei valley, as a result of the partial destruction of the dam and the lake outburst, a catastrophic mudflow arose, which traveled a path of approximately 10.5 km to the mouth of the river and was discharged into the Chuya river channel, the main water artery of the SE Altai. In the Maashei valley, the passage of mudflow led to the deposition of a huge amount of predominantly boulder-pebble material, to the destruction (in some places complete) of the larch forest in the floodplain, pruning of slopes, and in many areas to a change in the nature of the river channel. As our field observations in 2022 showed, the consequences of the mudflow that occurred on 15 July 2002 were clearly visible in the relief and sediments of the valley ten years after the event.

As a result of the outburst, the lake was completely drained (Figure 21).

5. Conclusions

Information about the current state of the glaciers of the North-Chuya range at different stages of its research since the first half of the 20th century, on the one hand, reflected the dynamics of glaciers, mainly their reduction, but on the other hand, reflected changes in research methods and an increase in the accuracy of the data obtained. Therefore, despite the general trend toward the glacier reduction after the maximum of the Little Ice Age, the number and total area of glaciers fluctuated with each cataloguing (

Table 10. Information about the glaciers of the North-Chuya ridge from different sources.

Number of Glaciers	Total Area, km2	Materials	Year	Source
118	127	Field observations	1920-s	[31]
116	136	Field observations	1940-s	[70]
168	143.16	Aerial imagery (1952, 1955), field observation of the 1960s	1960-s	[30]
-	154.09	Aerial imagery (1952) and field observations	1995	[73]
181	164.2	Data from the Glacial Inventory of the USSR, the rates of degradation in 1952–1995 extrapolated until 2003	2003	[33]
-	148.66	ASTER satellite imagery	2008	[22]
	112.9	Sentinel-2	2017	[34,74]
224	120.02	Sentinel-2	2021	This article

).

Most of the information in Table 10 is given for time points or periods that took place many years before 2021. Moreover, most of the data were obtained from field research, including interpretation of aerial photographs, in which distortions of the area and shape of glaciers are inevitable, or from ASTER images with a relatively low resolution. This is why the most justified is the comparison of our data with the results of the work on creating a unified catalog of Russian glaciers, where Sentinel-2 images were also used (https://www.glacru.ru/, accessed on 1 November 2022 [34,74]).

Significant differences between our results and the data of the catalog of glaciers in Russia are mainly related to the determination of the boundaries of glaciers in the areas of development of debris-covered ice and the assignment or non-attribution of the corresponding areas to the glaciers. In addition, over four years, some glaciers managed to retreat to a distance of up to 100 m (Figure 22).

The scale of glaciation during the LIA and its subsequent reduction were estimated for the North-Chuya range in the works of P.A. Okishev: the total area of the LIA glaciers was 192.18 km² (19.8% decrease by 1995) [73], i.e., 181 glaciers with the total area 208.2 (21% decrease by 2003). These estimates relate to a limited number of glaciers, mostly large ones, and do not affect the currently empty cirques, where glaciers have completely degraded with the LIA. This explains the large differences from our estimate: 400 glaciers with a total area of 304.90 ± 23.49 km² and a 61% decrease of the total area by 2021. The values of the average ELA changes (207 m) are also higher than the values given by Okishev for the largest glaciers of the North-Chuya ridge (65–70 m) [24].

The reduction in the area of glaciers of the North-Chuya range was somewhat larger than the average value obtained by us for 27 glacial centers in the southern part of Altai for the period from the ELA maximum until 2006–2021 (59%). The ELA rise was almost twice as high as that for the southern part of Altai (106 m). This fact corresponds to the regularity established by us for the southern part of Altai: in the more humid areas, the rise of the ELA was higher than that in the arid areas, which could be caused by a significant decrease in precipitation after the LIA maximum in the relatively humid parts of Altai, while in the more arid part, it changed little [26].

The results of our estimations of a decrease in the North-Chuya range glaciers after the LIA are quite comparable with the information from other regions of the world. In the European Alps, the total area of the glaciers decreased by 54.5% in the period of 1850–2000 [34]. Our North-Chuya ridge glaciers’ reconstruction for the same period had a 54% decrease. This similarity is quite understandable if we consider the position of both mountain systems at approximately the same latitude and the similarity of their altitudes. The estimates of ELA rise for the entire European Alps were also close in their value to what we obtained for the North-Chuya glaciers: ELA changes for the period between the mid-19th century and ~1970, based on AAR calculations, were evaluated ~69 m for the Swiss Alps [75] and ~94 m for the Austrian Alps [76], with a further rise of the ELA between 1984 and 2010 by about 170 m for the Western Alps [77], which gave the total ELA rise for the European Alps of over 200 m.

In other mountainous countries, the glacial shrinkage after the LIA maximum was less pronounced. In the Central and Western Himalayas, the total length and total area of 220 glaciers had decreased respectively by 35% and 31% by 2005–2015 since the LIA, the mean increase of ELAs from the LIA to the first decade of the 21st century, as reconstructed using the toe-to-ridge altitude method (TRAM), was 123 m [78]. In southeast Tibet, an average retreat of ~27% (length change) has been revealed, and there is a trend toward a stronger retreat for smaller glaciers, with an average rise in the ELA of ~136 m since the LIA [40]. It is likely that the less-pronounced changes in glaciers for these countries are associated with a greater average size of glaciers and their vertical range, which lead to their greater stability.

The main features of dynamics that we had revealed for the Bolshoi Maashei glacier in our research are the following: generally slow retreat of the glacier in 1850–1924, uneven rates of retreat between 1924 and 1962 with probable close-to-stabilization intervals in the late 1920s and in 1944–1947, a significant slowdown of retreat in 1962–1989, a faster retreat in 1989–2010, and rapid acceleration after 2010. The slowdown of the glacial retreat could be caused by the slightly positive mass balance trend as a result of lowering summer temperatures between 1920 and 1970. A rapid increase in summer temperatures in the 1990s changed the mass balance trend to negative: despite the stabilization of the summer temperatures’ level after 2000, the negative mass balance trend continued due to the decrease in the annual precipitation.

The acceleration of glacier retreat in the last one–two decades is in line with the global trend [16,79]. At the same time, regional climate changes have their own specifics, in particular, fluctuations in the amount of precipitation can differ significantly even for neighboring glaciation centers within the same mountainous country, introducing local features into the dynamics of glaciers.

The most well-studied glacier of the North-Chuya ridge and Altai as a whole, the Malyi Aktru, is located about 11 km to the south-east of the Maashei glacier. According to the multiyear mass balance observations [36] that were held between 1961 and 2012, there were years with a positive mass balance for the Maliy Aktru glacier, especially between 1965 and 1995. Observations of the glacial front variation there started in 1911 by Sapozhnikov, and in the 1950s they became regular. Unfortunately, after 2005 they were interrupted. The stabilization of the edge of the glacier took place around 1965–1970, and 1978–1981, and after 1991, the glacier accelerated its reduction, in some years to values of more than 40 m [80]. It is noteworthy that despite the general retreat of the glacier edge, until 1990, the volume of the glacier increased. It is likely that if it were not for the sharp warming in the 1990s, it would have gone on the advance. This circumstance agrees well with our mass balance calculations of a slightly positive mass balance trend between 1920 and 1970 for the Bolshoi Maashei glacier.

In other parts of Altai, the behavior of large glaciers has both similarities and differences with the dynamics of the Bolshoi Maashei glacier. The valley glaciers of the Tavan-Bogdo massif, which is situated about 100 km to the south, after the mid-1960s gradually accelerated their reduction, and a slight slowdown in retreat for a number of glaciers took place in the 1980s. In the 1990s, the glacier retreat accelerated, and the highest retreat rates are noted after 2010, reaching 40–70 m/year [26]. The acceleration of glacier retreat in recent decades was also caused by a sharp increase in summer temperatures in the 1990s, with further stabilization at high levels, noted by the Kosh-Agach, Altai (China), and Altai (Mongolia) meteorological stations [48]. However, weather station data did not indicate a decrease in precipitation in recent decades; therefore, after a sharp decrease in the mass balance of glaciers in the 1990s, after 2000, mass balance values have stabilized [81,82], which differs from the results of our calculations for the Bolshoi Maashei glacier, where it continues to decline.

For the valley glaciers of the Mongun Taiga massif as a whole, there was also a trend towards an acceleration of glacier retreat, especially pronounced from 2013 to 2016, when the two largest glaciers retreated by more than 40 m per year. Here, the warming of the 1990s was also the reason for the acceleration of the reduction, but it was combined with a decrease in the amount of precipitation (meteorological station Mugur-Aksy), so the decrease in the mass balance in these years was especially sharp. On the contrary, in the last 20 years, the amount of precipitation has been increasing and there has been some increase in the mass balance, which, nevertheless, has not yet allowed it to reach the level before the 1990s [82,83].

An abrupt decline of the mass balance in the 1990s followed by stabilization has been revealed for the Shapshal ridge [84], adjacent from the north to the Mongun-Taiga massif; this decline was caused by a simultaneous rise in the summer temperatures and a decrease in precipitation. However, no acceleration of the glacial retreat has been found for this glacial center, probably because there are no valley glaciers present there and cirque glacier Mushtuk, that has been studied, is less sensitive to climate changes.

No acceleration of glacier retreat has been revealed for the Tsambagarav massif glaciers (about 280 km to the south-east of the Maashei glacier) [82].

A review of regional data shows that, in general, the features of the dynamics of the Bolshoi Maashei glacier are characteristic of other large glaciers, in particular, this concerns the current sharp acceleration in the reduction of glaciers. In terms of mass balance, until the 1990s, the glacier was in more favorable condition in relation to the glaciers of the more arid regions of the eastern and south-eastern Altai. In the past 20 years, the climatic conditions for the existence of the glaciers of the North-Chuya ridge have continued to become less favorable, which should contribute to the high rate of retreat of its largest glaciers in the near future.

Here are some conclusions of our research:

Since the LIA, the glaciers of the North-Chuya range have been experiencing a reduction in area, while the relative reduction in the area of glaciers and the ELA elevation here exceeds the average values typical for the glacial centers of Altai.

A large reduction in the total area of glaciers before 2000 was caused by the predominant degradation of small glaciers on the reduced periphery of the ridge.

The retreat of the Bolshoi Maashei glacier after the LIA was slow most of the time. However, after 2010, the Bolshoi Maashei glacier dramatically accelerated its retreat, which is typical for most large valley glaciers in the region.

In the mid-1990s, there was an abrupt increase in summer temperatures, which greatly changed the mass balance of the glaciers of the North-Chuya range in a negative direction; in the last 20 years, the conditions for the existence of glaciers continue to deteriorate due to a decrease in precipitation.

The rapid degradation of glaciers reflects an increase in ablation and is accompanied by an increase in the risk of outbursts of high mountain lakes, which confirms the example of Lake Maashei.

The prospects for further work on this topic are seen as follows. It is necessary to obtain more complete information about the rate of reduction of the territory’s glaciers by obtaining information on additional time slices, for example, by using Corona images. It is required to resume field observations of the retreat of the largest glaciers of the ridge and restore the data series on their reduction by eliminating gaps using a combination of field and remote data. It is necessary to resume mass balance observations in previously well-studied areas (Aktru) and obtain mass balance information on representative glaciers of the ridge to predict the future behavior of glaciers.