Development of proglacial lakes and evaluation of related outburst susceptibility at the Adygine ice-debris complex, northern Tien Shan

The formation and development of glacial lakes in mountainous regions is one of the consequences of glacier recession. Such lakes may drain partially or completely when the stability of their dams is disturbed or as a consequence of impacts. We present a case study from the Central Asian mountain range of Tien Shan – a north-oriented tributary of the Adygine Valley, where the retreat of a polythermal glacier surrounded by permafrost has resulted in the formation of several generations of lakes. The aim of this study was to analyse the past development of different types of glacial lakes influenced by the same glacier, to project the site's future development, and to evaluate the outburst susceptibility of individual lakes with an outlook for expected future change. We addressed the problem using a combination of methods, namely bathymetric, geodetic and geophysical on-site surveys, satellite images and digital elevation model analysis, and modelling of glacier development. Based on this case of the glacial lakes being of varied age and type, we demonstrated the significance of glacier ice in lake development. Lake 3, which is in contact with the glacier terminus, has changed rapidly over the last decade, expanding both in area and depth and increasing its volume by more than 13 times (7800 to 106 000 m³). The hydrological connections and routing of glacier meltwater have proved to be an important factor as well, since most lakes in the region are drained by subsurface channels. As the site is at the boundary between continuous and discontinuous permafrost, the subsurface water flow is strongly governed by the distribution of non-frozen zones above, within, or beneath the perennially frozen ground. In the evaluation of lake outburst susceptibility, we have highlighted the importance of field data, which can provide crucial information on lake stability. In our case, an understanding of the hydrological system at the site, and its regime, helped to categorise Lake 2 as having low outburst susceptibility, while Lake 1 and Lake 3 were labelled as lakes with medium outburst susceptibility. Further development of the site will be driven mainly by rising air temperatures and increasingly negative glacier mass balance. All three climate model scenarios predicted a significant glacier areal decrease by 2050, specifically leaving 73.2 % (A1B), 62.3 % (A2), and 55.6 % (B1) of the extent of the glacier in 2012. The glacier retreat will be accompanied by changes in glacier runoff, with the first peak expected around 2020, and the formation of additional lakes.

1 Introduction

Glacial lakes often appear as gems gleaming amid harsh mountain environments. At the same time, they can pose a serious threat to downstream settlements and infrastructure. If we focus on Asian mountain ranges, the now common term “glacial lake outburst flood” (GLOF), first used for the Himalayan region, is ubiquitous. Mainly due to the development of satellite-based sensor technologies, the regions under the scrutiny of investigators have increased from the Himalayas (Benn et al., 2012; Fujita et al., 2009; Shrestha et al., 2010) and Karakoram (Chen et al., 2010; Haemmig et al., 2014) to the Tibetan Plateau (Liu et al., 2014; Wang et al., 2013; Zhang et al., 2017), Pamir (Mergili and Schneider, 2011), and Tien Shan (Bolch et al., 2011; Engel et al., 2012; Narama et al., 2017; Sorg et al., 2012). In the territory of Kyrgyzstan, there are about 2000 glacial lakes (>1000 m²), almost 20 % of which are potentially dangerous, and ∼15–20 lakes are at risk of sudden drainage each year (Erokhin and Zaginaev, 2016). However, this dry and glacier-meltwater-dependent region of Central Asia has remained rather in the background of research interest.

Some glacial lakes form and drain within a relatively short time (Erokhin et al., 2017). Others exist for years and decades without major change. Lastly, there are lakes that are in a phase of enlargement. The expansion of glacial lakes broadly correlates with climatic warming trends and negative glacier mass balance patterns (Lei et al., 2012; Mergili et al., 2013; Zhao et al., 2015). In addition, a recent study has shown that the expansion of glacial lakes under a dry continental climate regime can be closely related to thermokarst processes in permafrost as well (Li et al., 2014). This means that lakes enlarge when in contact with a retreating glacier tongue by the melting of ice in the lake basin bed or around the sides or by filling due to increased inflow from a glacier. Other factors that play a role in lake evolution include the erosion of surface drainage channels, the formation, expansion, or blockage of subsurface drainage channels, dam morphology changes, or slope movements adjacent to the lake.

Outbursts from mountain lakes occur almost annually in Kyrgyzstan (in recent years, these included Lake Merzbacher, 2017, Lake Chelektor, 2017, and Lake Aksai, 2015), and given the generally increasing use of mountain valleys, greater losses can be expected as a consequence (Dussaillant et al., 2010). Once an outburst is triggered, a flood (often evolving into a hyperconcentrated flow or debris flow) will threaten areas in the lower parts of a valley. To reduce the risk, an assessment framework (see, e.g. the GAPHAZ Technical Guidance Document; GAPHAZ, 2017) can be applied, which includes the following steps: (i) identification of potentially dangerous lakes, (ii) detailed evaluation of the hazard, and (iii) application of mitigating measures. Our study includes the second step – a detailed evaluation of the outburst susceptibility of selected lakes that had already been identified as potentially dangerous (Erokhin and Zaginaev, 2016). As the need for field data for such an assessment has been generally acknowledged, here, we present an assessment based on a combination of on-the-spot and remotely obtained data.

Although the study site is not particularly suited for the purpose of only monitoring glacier retreat, we did find it most appropriate in terms of observing the site in a more holistic way. The glacier itself is considered to be representative of smaller glaciers, which are subject to relatively rapid retreat rates (Narama et al., 2010) and generally have shorter response times to changes in climatic conditions compared to larger ones (Wang et al., 2014). This allowed us to follow a larger development period at a smaller timescale. Also, the existence of glacial lakes documenting glacier shrinkage was convenient – being at various stages of development, they were seen as able to help us to better understand the evolution of different types of lakes over time. Due to the occurrence of lake outbursts in the region, the high elevation of the site, and its proximity to populated areas, an outburst susceptibility assessment strategy, with an estimation of future change, was considered to be desirable and of practical relevance.

Based on a detailed study of one particular site, we looked at the issue of glacier development and the associated, actively evolving glacier forefield and its consequences on downstream areas. We have used the term “ice-debris complex” (as coined by Bolch et al., 2019) to describe the glacier and the debris landforms with varying shares of buried ice in front of it. In order to address the problems associated with the rapid processes that operate in mountainous regions, we incorporated perspectives of future scenarios for the development of the site, and the related outburst susceptibility, into our study. Because our intention was to promote a more complex approach to the issue of susceptibility assessment, our study included three particular objectives: (i) to analyse the past evolution of different types of glacial lakes influenced by the same glacier, (ii) to provide insight into the probable future evolution of the site, and (iii) to evaluate the outburst susceptibility of individual lakes, with an outlook for possible future change. As a consequence of our objectives, the structure of the paper is as follows: first, we present an analysis of the past development of the local setting, following both the glacier and the proglacial lakes; then, we report on the potential future conditions for the site; and, finally, we include evaluations of the lakes' outburst susceptibilities and estimations of how they will likely evolve under changing climatic conditions. In our opinion, to analyse glacier retreat and lake formation without touching on the topic of consequences (present and near future) would be as incomplete as evaluating outburst susceptibility without knowing the past evolution of the site.

2 Study area, data, and methods

The research site is situated in northern Tien Shan, specifically Kyrgyz Ala-Too. The highest peak of this east–west-oriented range is Pik Semenova Tienshanskogo, at 4895 m a.s.l. Although the Kyrgyz Ala-Too is not as severely glaciated as the ranges in central Tien Shan, the currently glacier-free parts of the valleys do show signs of the former glacier extent. The main Ala Archa valley currently has ∼33 km² of its area covered by glaciers. It drains into the Chuy River, which is part of a large endorheic basin. The most common lake type of Kyrgyz Ala-Too are lakes formed in intramorainic depressions (Table 1). According to Erokhin (2012), 83.6 % of potentially dangerous lakes are ice-cored moraine-dammed lakes, which include intramorainic and thermokarst lakes in the presence of buried ice.

Table 1 Representation of different lake types in Kyrgyz Ala-Too, northern Tien Shan, in 2017. Only lakes with a minimum area of 1500 m² were categorised. Based on manual mapping using satellite imagery in Google Earth and refined with the field data included in Erokhin and Zaginaev (2016).

The development of glaciation in the Ala Archa Basin is monitored, in the long term, mainly due to its position in proximity to the Kyrgyz capital, Bishkek, and the popularity of the valley with tourists visiting the National Park. According to Aizen et al. (2006), there was a change in the glaciated area of the Ala Archa watershed by −0.12 km² yr⁻¹ (−0.29 % yr⁻¹) between 1963 and 1981, and by −0.2 km² yr⁻¹ (−0.51 % yr⁻¹) between 1981 and 2003, documenting an accelerated retreat by the end of the 20th century. The overall reduction in glaciated area over 60 years (1943–2003) was 15.22 % (a loss of 6.52 km²; Aizen et al., 2006), which is an above-average value; in Tien Shan, the glacier area shrank by 14.2 % over the same interval (Aizen et al., 2007a). Farinotti et al. (2015) confirmed this trend, and Bolch (2015) pointed out that, although the number of glaciers increased due to the disintegration of several larger glaciers, the total glaciated area reduced by 18.3±5.0 % between 1964 and 2010. Rising air temperatures (especially since the 1970s) have caused a negative mass balance in most glaciers, which has also led to changes in the hydrological regime of the glacial streams (Aizen et al., 1996; Glazirin, 1996; Pieczonka et al., 2013). Sporadic permafrost occurs at elevations above 2700 m a.s.l., is discontinuous at 3200–3500 m a.s.l., and is continuous above ∼3500 m a.s.l. (Gorbunov et al., 1996).

The Adygine watershed has an area of 39.6 km², and its glaciated area constitutes ∼10 % of the basin (3.9 km²). The total elevation differs by up to 2370 m; the valley below the glacier has an average slope of 10.8^∘. The Adygine ice-debris complex (42^∘30^′10^′′ N, 74^∘26^′20^′′ E) closes this 8 km long tributary valley with a northern orientation and reaches an elevation of 3400–4200 m a.s.l. (Fig. 1), which means it lies mainly within the zone of continual permafrost (Gorbunov et al., 1996). The permafrost is known to reach depths of up to 100 m, and the active layer can be up to 4 m thick at similar altitudes in northern Tien Shan (Marchenko et al., 2007). The occurrence of ice-rich, perennially frozen ground thus controls the underground water flow at the site, allowing it only to pass through non-frozen zones (mainly the active layer and taliks; Cheng and Jin, 2013). The large landforms of mixed ice debris below the glacier were formed by long-term viscous creep (Bolch et al., 2019) and contain buried glacier remnants and ice lenses. Several such ice lenses were recently exposed due to surface erosion and may have had an impact on lake evolution and stability (Karlsson et al., 2012).

Figure 1 Overview of the Adygine ice-debris complex (a) and its recent dynamics; exposure of buried ice embedded within perennially frozen ground (b, c). 1–7 are the studied lakes; P1–P3 refer to geophysical profiles. Study site marked with an arrow on map of Kyrgyzstan (upper left corner). Positions of (b) and (c) are marked on Fig. 1a. Image source: Worldview-2 (2011); photo: K. Falatkova (2017).

The upper part of the complex involves a polythermal glacier (2.8 km²), followed by a three-level cascade of glacial lakes that have evolved as a consequence of glacier retreat over the past 50 years. The glacier terminus is currently situated at 3600 m a.s.l.; the tongue is rather steep and short, emerging from a relatively flat, larger source area. A large part of the glacier is covered with a thin layer of fine material, considerably lowering the albedo of its surface. The position of the equilibrium line altitude (ELA) on the northern slopes of Kyrgyz Ala-Too is at ∼3900 m a.s.l., resulting in an ablation zone covering more than 65 % of the glacier area. The lakes found at the site (Table 2) are of varying ages and have different positions in the connected hydrological system that drains meltwater down the valley. The most recently formed lakes are close to the glacier terminus, reacting sensitively to changes in glacier melt rates and redistributing the water further downstream, either by surface or subsurface channels. The middle level is represented by the largest lake of the site – Lake 2 (32 000 m²) – that has a fairly stable annual hydrological regime (Falatkova et al., 2014), supplying water to the lowest part of the cascade – Lake 1. Not being a permanent lake, this intramorainic depression is only filled with water during an ablation season, when the rate of incoming meltwater is higher than the capacity of its subsurface drainage.

Table 2 Basic parameters of the three largest lakes of the site.

¹ Potential volume of the lake in the case of drainage blocking.
² During observed high water level dam overflow, water sinks underground after approx. 5 m.

A Holocene moraine complex, at an elevation of 3450 m a.s.l., forms the lower part of the Adygine complex. It consists of several parts of varying age (Shatravin, 2000), some of them being identifiable based on visual inspection only. The western part is composed of creeping ice-rich debris, with a typical relief of arcuate ridges and furrows in its terminal part; ice lenses are present and occasionally exposed (Fig. 1b). From the east, a smaller ice-debris matrix moves under the influence of gravity and adjoins the central landform. The main body is a so-called debris-derived rock glacier that formed below the glacier terminus. Lying probably within the discontinuous permafrost zone, this landform has a prominently oversteepened front that is almost 700 m wide. In its western part, thermokarst processes are manifested in subsidence craters, cracks, and sagging, and numerous thermokarst lakes and exposed ice lenses (Fig. 1c) can be found here. The eastern part is considered to be of an older generation, with a rather flattened surface and stable lakes with lower turbidity.

Climatic conditions in the area are continental, characterised by relatively low precipitation and high annual and daily air temperature fluctuations. This is clearly evident at the Ala Archa station (2200 m a.s.l.), which reports a 22 ^∘C difference between the mean air temperature of the warmest and coldest months (13 ^∘C for July, −9 ^∘C for January), a mean annual air temperature (MAAT) of 3.3 ^∘C, and precipitation totals of 450 mm. Some 1300 m higher, at the Adygine site, the MAAT is around −3.5 ^∘C (based on weather station data from 2008 to 2013), and annual precipitation is estimated to be 700–800 mm. On the northern slope of Kyrgyz Ala-Too, the winters are characterised by low precipitation totals, whereas the maximum precipitation is typically reached in late spring and the beginning of summer.

2.1 Field mapping

A research station, built at the site in 2008 next to the Lake 2 surface outflow, provided an adequate location for field measurements to be taken. Climatic data and ground temperatures were obtained from an automated weather station situated near the Adygine glacier terminus at 3560 m a.s.l. The hydrological regime of the lakes was monitored with pressure sensors installed in Lake 1 (2012–2015), Lake 2 (2007–2017), and Lake 3 (2012–2017).

The glacier retreat has been monitored since 2006, during the ablation season, by means of a geodetic survey carried out with a Leica TCR 705 total station (infrared beam on a reflective prism, accuracy: 0.005 m). The survey was focused solely on the terminal part of the tongue, where the retreat is most significant and can result in the exposure of overdeepenings (potential future lakes). This time series was supplemented with glacier limits acquired from aerial images taken in 1988 and 1962 (scale of the survey: 1:38 600, image resolution: 1 m). The lakes' shoreline changes were measured by the same means as the glacier terminus. To capture changes in the lake basins and derive the volume of retained water, a repeated bathymetric survey was carried out at selected lakes (Sobr and Jansky, 2016). Depth was recorded in 2.5 m steps along defined profiles, using an echo sounder mounted on a boat (depth measurement accuracy: 0.1 m). We processed the data obtained and interpolated them into a bathymetric map.

The selected key locations at the site (the lakes' surroundings and their dams) were surveyed by electrical resistivity tomography and spontaneous polarisation, in order to detect buried glacier remnants, frozen ground, and possible seepage routes. In order to acquire information on glacier thickness and bed topography for adjusting the model of future glacier extent and identifying potential areas of lake formation, ground-penetrating radar (GPR) was used. The location of the profiles is shown in Fig. 1a.

The method of electrical resistivity tomography (ERT) consists of measuring the resistivity of the subsurface by means of a number of electrodes located along a profile. Interpretation of such measured data was performed using the software Res2Dinv (Loke and Barker, 1995), which provided 2-D inversion of the measured resistance data and calculated the 2-D cross section under the profile. The resistive section of the observed profile approximates the actual distribution of the individual resistive layers in the depth range given by the maximum spacing of the electrodes. We used the ARS-200E model (GF Instruments, CR), with 48 electrodes at 4.4 m increments. The maximum range of the electrodes was 206.8 m, with a depth range of about 50 m.

The spontaneous polarisation (SP) method measures the natural electrical potential of the rock environment; in the case of lake dams, the filtering potential that arises from water filtration through the porous environment can be tested. Measurements were made with a GEOTOR I with offset compensation. Non-polarisable electrodes were used to minimise the effect of transient resistance between the measuring electrode and the geological environment. The measurements were carried out using potential variation, in which the electrical potential of a given point on the profile is measured relative to a fixed reference point. The measurement steps were around 10 m, and the location of every fifth point was recorded using GPS.

GPR was used to acquire information on the glacier thickness and subglacial topography. GPR data were collected along the central longitudinal profile, for five transverse and three short, connecting profiles. An unshielded 50 MHz rough terrain antenna and MALÅ CU-II control unit were used for data collection. The signal acquisition time was set to 860 ns and the scan spacing was about 0.07 m. The collected data were processed and interpreted using Reflex-Win software version 4.5.

2.2 Modelling of glacier evolution

The glacier mass balance model – glacier evolution runoff model, or GERM – applied in this study has been well established for simulating glacier mass balance from climate data, based on a sophisticated degree-day approach (Hock, 1999; Huss et al., 2008). In our case, the future glacier mass balances of the Adygine glacier until 2050 were calculated with the GERM, forced by empirically downscaled scenarios' data for daily air temperature and precipitation, as described below. As the glacier mass balance was not measured at the Adygine glacier, we used the mass balance data from the nearby Golubina glacier (WGMS, 2017) for validating the simulated mass balances of the GERM against measured mass balances for the period 1981–1994. Mass balances were validated for both winter and summer. In addition, we were able to use discharge data from the Adygine River (from gauging near the confluence with the Ala Archa River) for 1960–1987 to validate the modelled discharge of the GERM with observed data.

Despite the generally poor availability of data in Central Asia, daily air temperature means and precipitation totals were gathered from two stations close to the study area: Alplager (2130 m a.s.l., 1978–2010) and Baityk (1580 m a.s.l., 1914–1979). The meteorological station data were used to calibrate and validate the downscaling model (multivariate regression and analogue method; Benestad, 2004; Benestad et al., 2007). Reanalysis data from NCEP/NCAR (Kalnay et al., 1996) were used as predictor fields in order to downscale long-term datasets of precipitation and air temperature for the study region. For scenario simulations until 2050, three different SRES scenarios – A1B, B1, and A2 (IPCC, 2007) – were performed by the ECHAM5/MPI-OM coupled atmosphere–ocean model (Röckner et al., 2003). Although the SRES scenario data were replaced by the Representative Concentration Pathways approach of the Intergovernmental Panel on Climate Change (IPCC, 2013), the uncertainty ranges of future large-scale temperature and precipitation changes provided by the two approaches are rather similar and smaller compared to the uncertainty range of climate data coming from further downscaling. We argue that, for the purpose of assessing future glacier extent, in particular for the future development of glacial lakes, these scenarios are sufficiently realistic.

In order to estimate the potential locations for glacial lake formation, the data on glacier bed topography acquired by GPR were used (see Supplement). We also consulted the digital elevation model (DEM) data for the glacier surface from the Shuttle Radar Topography Mission, which has a resolution of 30 m (Farr et al., 2007). The combined information on future glacier extent and the topography of the exposed area was analysed in GIS by determining flow directions and identifying sinks (i.e. overdeepenings).

2.3 Outburst susceptibility assessment

The presented lake outburst susceptibility assessment is qualitative and was adapted to accommodate the distinctive regional features – dry climate, the common occurrence of buried ice and permafrost. We based the assessment procedure on knowledge we gained through fieldwork in the region, focused especially on the hydrological systems of proglacial lakes. The specificity of the region is the origin of the lakes – the most common type is lakes formed in intramorainic depressions (Table 1). Many lake hazard evaluations are based on cases of moraine dams with a distinct morphometry; however, this type of lake does not occur that frequently in the region, representing only 2 % of the potentially dangerous ones (Erokhin and Zaginaev, 2016). Nevertheless, inspiration was drawn from assessments by ICIMOD (ICIMOD, 2011), Allen et al. (2016), Frey et al. (2010), and Huggel et al. (2004). Hazard assessments based solely on remotely sensed data are appropriate when applied over a large area and as a first step; however, at least for lakes labelled as potentially dangerous, additional evaluations, which would include on-site investigation, need to be carried out.

In this paper, the overall susceptibility is introduced as being a combination of conditioning factors (i.e. a lake's inherent characteristics determining the lake's susceptibility to bursting and causing flooding) and the presence of possible triggers that would have the capacity to cause an outburst. To assess a lake's inherent susceptibility, the following parameters were selected and a simple qualitative rating suggested (Table 3):

• Lake volume involves the size of the lake and its volume of retained water. Although outburst susceptibility does not increase proportionally with lake size, a larger volume generally means greater hydrostatic pressure on the dam and an increased potential to cause damage. The thresholds were set according to information on GLOF cases in Kyrgyzstan.

• Lake type is primarily connected to the stability of the dam. The material forming the dam (or the depression in which a lake exists) has varying characteristics and behaviours in relation to water. The type of rock step is considered to be the most stable, as water usually only has a minor effect on this; however, when covered with the loose morainic material left behind a retreating glacier, this part of the dam can be prone to erosion and piping and will become less stable overall. The interaction of a lake with an ice dam results in even higher risks – ice melting due to heat transfer from water and high hydrostatic pressure on the ice dam leading to its uplift, causing lake drainage.

• Ice contact is another parameter that influences lake stability. A lake in direct contact with a glacier is considered to be more dangerous than one at a distance, as it will be strongly affected by glacier behaviour. Its evolution and morphological changes can be dynamic as the glacier retreats, advances, or disintegrates.

• Drainage type is the parameter that represents a lake's hydrological regime, covering only two options – surface or subsurface drainage. It is the limited (or changing) capacity of a subsurface channel that makes lakes with this type of drainage more dangerous, as such lakes are prone to filling when inflow is increased. Accumulated water can then be released rapidly due to channel enlargement under higher hydrostatic pressure. A surface channel, in comparison, regulates the lake water level naturally and maintains its hydrological balance.

• The growth possibility of a lake depends on lake basin characteristics and the presence of ice that can melt. Lakes with terminated expansion possibilities are considered more stable than those which can enlarge for some reason.

Table 3 Inherent susceptibility parameters and ratings of lakes. The scheme on the right shows possible combinations of parameter ratings (H – high, M – medium, L – low) and the resulting susceptibility. Note that this simplifying table serves as a guideline for the susceptibility estimate and that even medium resulting susceptibility means a lake is dangerous and may burst and cause damage.

The determination of the degree to which a lake is susceptible to outburst is followed by an assessment of potential triggers. The means of determining whether a trigger has the potential to cause outburst in the specific lake or not may be distant (satellite imagery, DEM), as is the case for the first group of triggers in Table 4, or can be on-site, as is necessary for inner dam stability or subsurface channel functioning, for example.

Table 4 Triggers with potential to cause lake outburst and their description.

¹ Avalanche starting zone steepness: polythermal glacier over 25^∘ and cold-based glacier over 45^∘ (Alean, 1985).
² Average slope trajectory over 14^∘ (Noetzli et al., 2003).

The fall of mass into a lake is one of the most common triggers of lake outburst (Ding and Liu, 1992; Emmer and Cochachin, 2013; Falatkova, 2016), resulting in an impact wave that can overtop or destabilise the dam. Following a widely used procedure, the distance between a lake and a slope and slope steepness are the parameters that determine the trigger potential (Alean, 1985; Fischer et al., 2012; Noetzli et al., 2003). In the case of calving, a lake's contact with a cliff-like glacier terminus is a crucial condition.

The development and stability of a glacial lake can be affected by the melting of buried ice within a dam. During stagnation of a debris-covered glacier, an ice-cored moraine dam can be formed. These cores present a weak point in a dam, as their ablation changes the dam's inner structure (Richardson and Reynolds, 2000). Besides that, dam stability is also a function of the overall frozen or unfrozen condition of the dam material.

Not only landslide dams are prone to failure caused by internal erosion (piping). Cases of GLOF due to piping in moraine dams have also been observed (Xu, 1988). Moraines consist of unconsolidated heterogeneous material and therefore may be prone to seepage and piping as a conduit grows and a dam is weakened (Awal et al., 2011). Another dam-stability-decreasing phenomenon – earthquakes – has been confirmed as being a primary trigger of several lake outbursts (Clague and Evans, 2000; Lliboutry et al., 1977).

In the case of a lake not having surface drainage, and meltwater from a glacier being routed into the lake, there is the potential for a significant hydrostatic pressure increase at times of higher melt rates of snow or ice. A sudden air temperature rise could therefore cause a disturbance in the lake's hydrological balance, with possible consequences for dam stability. An increase in hydrostatic pressure can also be caused by the blockage of a subsurface channel, which is a very unpredictable, and still rather poorly understood, phenomenon. Monitoring changes in lake morphometry due to frequent sliding of lake walls or disruptions in the hydrological regime of a lake can serve as a tool for determining channel blockage potential.

Changes in lake outburst susceptibility result from a complex process associated with the interaction of geomorphological processes and the behaviour of the glacio-hydrological system. Development of the englacial and subsurface routing of meltwater is especially crucial; while being difficult to determine, this was beyond the scope of this paper. The presented estimate of susceptibility change is linked to climatic changes (rising MAAT) only and involves impacts accompanying deglaciation that have been observed in other mountain ranges, such as glacier retreat, altered glacier melt rate, slope instability (debuttressing effect), permafrost degradation, and buried ice melting (Allen et al., 2016; Frey et al., 2010; Haeberli et al., 2017). The general assumptions for future climatic conditions, glacier retreat, and altered melt rates are supported by the glacier mass balance model (Sect. 3.3).

3 Results

3.1 Past glacier development

The terminus position of the investigated glacier has changed significantly since the 1960s. Over the past five decades, the furthermost part of the glacier terminus has retreated ∼630 m (Fig. 2). The Adygine glacier probably became disconnected from its smaller western branch during the 1950s. According to an aerial image taken in 1962, the tongue was divided into two parts by a resistant rock outcrop. The larger western part reached the rock riegel at an elevation of ∼3540 m a.s.l. (close to the current Lake 2 drainage channel); the narrower eastern tongue bypassed the outcrop, reaching 100 m lower to 3450 m a.s.l. (currently the position of Lake 1). In the following decades, both parts retreated at an average rate of 8 m a⁻¹ in 1962–1988 and 14 m a⁻¹ in 1988–2006. Currently, there is no elevation difference between the western and eastern parts of the terminus (3640 m a.s.l.). In the latest period of observation (2007–2017), the terminus retreat rate had similar values to the earlier (1962–1988) phase, at slightly over 8 m a⁻¹. A significant recession period in recent years was recorded between 2013 and 2015, when the terminus retreated some 10–15 m a⁻¹. The positional differences of retreat rate are linked to the variability of the glacier bed topography (rock outcrops) and to heat transfer with meltwater (contact with glacial lake).

Figure 2 Retreat of the Adygine glacier tongue between 1962 and 2017. 1–7 are the proglacial lakes. The arrow marks the rock outcrop (riegel). Base image: Worldview-2 (2011).

3.2 Formation and development of the lakes

There are numerous lakes of varying size, type, and stage of development in the study area. They form a three-level cascade, and at least some of the lakes are hydrologically connected.

3.2.1 Lake 1

The lowest part of the cascade is formed from about a dozen thermokarst lakes, situated on a large morainic, perennially frozen platform at an elevation of ∼3450 m. The largest lake (Lake 1; Janský et al., 2010) is situated directly below a rock riegel, in a depression covering ∼200 000 m³, which was revealed when the eastern part of the tongue receded. The exact year when the depression filled with water is unclear. According to our knowledge of the drainage system, it was after Lake 2 developed and started supplying Lake 1 with a larger volume of meltwater (probably during the late 1980s). This intramorainic lake is not permanent, as the inflow only exceeds the outflow channel capacity during an ablation season. The lake level fluctuates significantly during the course of a day (leading to volume changes of between 10 000 and 15 000 m³), caused by the limited capacity of the subsurface outflow channels. Besides the surface water intake, water is probably routed there by seepage from Lake 2 and Lake 3. There have been changes in the morphometry of the basin in recent years, as sliding of unconsolidated material from the basin walls has been observed. The part close to the lake drainage has changed recently, extending the basin in a north-easterly direction.

3.2.2 Lake 2

Lake 2, the largest lake at the study site, began forming in about 1960, after the western glacier tongue retreated behind a rock riegel at 3540 m a.s.l. Since then, the lake area has extended, together with further terminus recession, reaching 32 700 m² in 2005 (first on-site survey). In the second half of the 1990s, the lake lost direct contact with the glacier, and, as a result, the lake ceased to grow. The lake volume is regulated by a surface outflow channel, eroded in morainic material, covering a rock outcrop. In the last decade, a slight change in the lake's area has been recorded, caused by siltation near the inflow. In 2017, the lake area was 30 900 m², and the maximal lake depth has also slightly decreased since the start of monitoring (from 22.2 m in 2008 to 21.3 m in 2015), probably also as the result of siltation. The lake volume changed accordingly from 208 000 m³ (2008) to about 195 000 m³ (2015). Geophysical sounding confirmed the presence of buried ice in the western part of the lake dam area – the remnants of a glacier tongue formerly bypassing the rock outcrop, now covered by 8–10 m of ice-rich frozen morainic material (Fig. 3). The overall thickness of ice and frozen ground exceeded the instrument depth range, being more than 40 m. The low resistivity values (<5000 Ω⋅m) marked a thin, active layer at the surface but also non-frozen material below the lake surface outflow and possibly the bedrock. Seepage was detected at the outflow (through the non-frozen ground above the bedrock) and at the western end of the dam, which can be linked to the active layer and non-frozen zones within the permafrost.

Figure 3 Electrical resistivity results along Profile 1 and Profile 2 for Lake 2 and Profile 3 for Lake 3. Resistivity values of materials depend on grain size, pore size, water content, salinity, temperature, and phase resistivity increases due to the change from electrically conductive water to non-conductive ice (Kneisel et al., 2008). For location of the profiles (P1–P3), see Fig. 1a.

3.2.3 Lake 3 and newly emerging lakes (4–7)

Since 2004, several lakes have formed in the proximity of the retreating glacier terminus, filling rather shallow depressions. The lakes' development is closely linked to the amount and routing of glacier meltwater as well as to basin stability. Our data suggest the importance of direct contact with the terminus, as only those lakes showed an increase in area. Lake 4 is an example of one that was enlarging slightly until it lost contact with the glacier in 2011–2012. Most lakes formed, enlarged, or remained stable over a period of several years and finally drained by subsurface routes (for detailed information on the lakes' morphometric changes, see the Supplement). An exception is Lake 6, which seemed stable over the monitoring period, until it suddenly drained in 2015 but was filled again the following year. The area exposed following the glacier retreat is underlain by buried ice lenses (visual inspection; Fig. 1a) within frozen ground (documented by the measured resistivity; Fig. 3). A possible explanation for the lake's sudden drainage might be permafrost degradation caused by warming by the lake water (Shur and Jorgensen, 2007).

Although these proglacial lakes are highly unstable, the hazard they pose is negligible, as only a small amount of water (up to 4000 m³ in Lake 4 and Lake 6) can be retained in such shallow depressions. Lake 3 is an exception; this easternmost lake is in direct contact with the glacier and grows annually. It is situated in a deep basin – in 2007, a bathymetric survey recorded a maximum depth of 3.7 m; in 2017, as the lake enlarged, it was over 14 m deep. Besides the climate-driven glacier retreat, the lake has been growing due to increased glacier melting, caused by heat transfer between the lake water and the ice. The lake expanded from 5710 m² in 2007 to 8830 m² in 2012. In the past few years, the growth has accelerated, leading to a lake area of 16 020 m² in 2017 (Fig. 4). Shortly after its formation, the lake had a surface outflow; however, the drainage soon switched to subsurface channels. Currently, during times of high lake water levels, water overflows the dam. This minor surface drainage sinks below the surface after 5–6 m, however. The geophysical survey (Fig. 3) revealed the dam structure, which is mainly composed of frozen debris with ice lenses. In the eastern part, there is distinct glacier ice covered with ∼10 m of ice-rich frozen material. About 15–20 m below the lake outflow, there is a zone of lower resistivity (several thousand ohm-metres, Ω⋅m), which could indicate either bedrock or a non-frozen zone (talik). If the latter is the case, the lake's further deepening, and expansion of the non-frozen zone beneath the basin, could result in a connection being formed with this talik. As the volume of the retained water is increasing rapidly (7800 m³ in 2007, 29 300 m³ in 2012, and 106 000 m³ in 2017), this recently insignificant lake has become a focus of attention.

Figure 4 Bathymetric map and area changes of Lake 3 between 2007 and 2017. Lake depth data obtained during a field survey at the site carried out in 2017.

At the level of Lake 1, there are several small glacial lakes of varying age (Fig. 5), which were formed by thermokarst processes. As these lakes only retain small volumes of water, and no significant areal development has been observed, a closer investigation from the outburst susceptibility point of view was not carried out. However, a few lakes situated in the youngest generation of a moraine recently uncovered the basin sides formed by buried glacier ice (Fig. 5, left). These lakes are expected to develop further in the future, either by enlarging or draining through newly opened subsurface channels.

Figure 5 Thermokarst lakes of varying age in the moraine complex below the Adygine glacier. These small lakes are situated at ∼3450 m a.s.l., all within a few hundred metres of each other. Differences in lake water turbidity or colour suggest existing (left) or missing (right) recharge from melting buried ice. Photo: K. Falatkova (2017).

3.3 Expected future site conditions

According to all three modelled scenarios, an increase in air temperature in the area is highly likely in future decades, with an estimated higher temperature increase in spring and summer and a moderate temperature increase in autumn and winter. Scenarios of annual precipitation showed a reduction under all three scenarios. In particular, summer and autumn depicted a decrease in precipitation. Winter and spring displayed a minor reduction, or even a slight increase, in precipitation totals (for winter under A2, for spring under B1). Consequently, future mass balances simulated for the Adygine glacier by the GERM are negative, and runoff will be altered. Low-pass-filtered values of the mass balance of the Adygine glacier were negative throughout the simulation period (Fig. 6a). Scenario A1B showed significantly less negative mass balances between 2015 and 2035 compared to the other scenarios. In this period, Scenario A2 was characterised by the highest decadal variation, with filtered values of between −1300 and −500 mm w.e. After 2030, all scenarios illustrated a significant negative trend, indicating that a tipping point for glacier shrinkage will be reached. Scenario A1B showed the steepest decline for this period, with a gradual recovery after 2045.

Figure 6 Glacier mass balance (a) and glacier areal evolution (b) modelled according to Scenario A1B, Scenario A2, and Scenario B1 for the period 2010–2050. Data smoothed using 10-year low-pass filter.

Due to the negative mass balances (Fig. 6a), all scenarios resulted in a reduction in glacier area. Under Scenario A1B, the glacier area reduced only slightly in the first two decades (2015–2035), when compared to the other two scenarios (Fig. 6b). After 2035, Scenario A1B showed accelerated areal shrinkage (resulting from the accelerated mass loss). All three scenarios simulated a significant decrease in glacier area by 2050, specifically leaving only 73.2 % (A1B), 62.3 % (A2), and 55.6 % (B1) of the 2012 glacier area. The model simulations indicate (Fig. 7) that, under Scenario B1, the glacier will disintegrate into several parts, whereas, under Scenario A2, only one smaller part will remain. The modelled glacier retreat implies the potential for new lakes to form in the exposed area. The glacier bed topography (Fig. A2), based on GPR profiling, revealed several overdeepenings. Under Scenario A1B, three new lakes, with a total area of ∼0.1 km², could emerge, while Scenario A2 presented the potential for seven new lakes to form, with a total area of 0.13 km². Scenario B1 exhibited the greatest modelled glacier retreat, and the exposed terrain had the potential for 11 new lakes, with a total area of 0.17 km².

Figure 7 Modelled future glacier extent and potential areas for lake formation for 2050 under Scenario A1B, Scenario A2, and Scenario B1.

As a result of greater negative glacier mass balances, total runoff from the Adygine glacier catchment is expected to rise. According to Scenario A1B and Scenario A2, the total runoff will increase by 19.7 % and 25 %, respectively, between 2010 and 2050 (Fig. 8a). However, Scenario B1 showed a rather sensitive reaction of the glacier dynamics to climate change, with a clear peak in discharge values around 2020. This was followed by a decline and overall stagnation in the following decade, resulting in no significant trend in the modelled period.

Figure 8 (a) Accumulated annual runoff and monthly runoff for May between 2010 and 2050. Data smoothed with a 10-year low-pass filter. (b) Mean daily runoff (averaged over all scenarios) from snowmelt and ice melt for the periods of 2010–2020 and 2040–2050.

Obviously, runoff for individual months will change in the future. The most pronounced changes were seen at the beginning and end of an ablation period – specifically, in April, May, and October. Slight increases in runoff are also expected in June and September. However, for the core period of the ablation season – July and August – no significant changes were modelled. Separating total runoff into its two main components – melt from snow and melt from ice – for both today and the future (Fig. 8b) highlighted a significant temporal shift for the snow component. The time of peak runoff from snowmelt, which was in mid-June for the first decade of 2010–2020, will change to mid-May in the decade 2040–2050. This will coincide with a slightly earlier start of ice melt at the glacier, as the snow cover will have depleted earlier as well, thus exposing bare ice.

Table 5 The resulting lakes' inherent susceptibility to outburst.

3.4 Outburst susceptibility assessment and its expected change

The inherent outburst susceptibility of the studied lakes, based on the lakes' characteristics, determined here by five parameters, is summarised in Table 5. Lake 1 and Lake 3 were marked as having medium susceptibility, whereas the largest lake of the site, Lake 2, was considered to be less inclined to burst and cause a flood. Several areas were determined as being a potential source of mass movement that could reach the lakes (Fig. 9). These are mainly the steep upper part of the glacier with visible deformations and cracks and also the eastern lateral moraine with buried ice. Both surface and subsurface hydrological connections among the lakes indicate routing of excess water in the case of rapid melting or lake drainage. As shown by the model projections, and also indicated by the position of the ELA, the glacier is expected to shrink significantly in the coming decades. Sorg et al. (2014) confirmed that future susceptibility is likely to change due to glacier retreat and the accompanying consequences. For instance, the slopes may become more prone to destabilisation due to permafrost degradation, in particular those that are south-oriented (i.e. exposed to direct solar radiation).

Figure 9 Geomorphological conditions and potential triggers of lake outburst at the Adygine ice-debris complex. View from the north (Google Earth, 2018).

Lake 1 has a medium susceptibility to cause flooding, mainly due to its origin (intramorainic depression), subsurface drainage, and the potential for enlargement. Outburst could be triggered by a fall of mass from the north-oriented slope of unconsolidated debris, which contains buried ice. The weak point is the subsurface drainage channel that leads through the moraine complex. Changes in the capacity and, above all, the blockage of such a channel would lead to the lake filling up and subsequent outburst. That could also be initiated by the outburst or overflow of Lake 3. Its outburst susceptibility is expected to increase in the future, due to the intensified melting of buried and exposed ice near the lake basin and destabilised slopes in its vicinity.

Lake 2 has a low susceptibility to cause flooding because of its stable surface drainage, distance from the glacier terminus and the associated ending of lake basin development. As there are no proximal unstable slopes (only rock outcrops), outburst is not likely to be triggered by a displacement wave from the fall of mass into the lake. A possible weak point is represented by buried ice and water seepage in the western part of the dam that could lead to changes in the subsurface drainage network. A change in outburst susceptibility is expected to be negligible in the future, as surface drainage would be able to deal with possible greater inflow, and the lake stability will not be threatened by unstable slopes.

Lake 3 has a medium susceptibility to cause flooding, due to it being in contact with the glacier, having relatively unstable subsurface drainage, and having the potential to grow both in area and depth. Despite the contact with the glacier tongue, a dam overflow by impact wave from calving is not probable because the terminus is not crevassed and a calving front has not developed. However, the lake is within reach of an ice avalanche from the steep glacier slope and also a landslide from the adjacent, unstable lateral moraine. The eastern part of the dam is prone to failure, due to the presence of buried ice. Overflow of the dam caused by increased inflow from the glacier may lead to progressive dam erosion. Outburst susceptibility is expected to increase in combination with increasing lake volume and hydrostatic pressure on the dam. Also, destabilised slopes in the lake's proximity will threaten its stability.

Individually, the lakes present a certain amount of threat. Nevertheless, as these lakes are interconnected, the overall hazard represented by the site should be considered. In this case, a chain reaction is a very probable scenario. The lowest lake (Lake 1) is subject to changes in its basin morphometry and thus in the functioning of its subsurface drainage channel. In the case of a rapid increase in drainage from Lake 2, or outburst from Lake 3, this lake would be hit, although it is rather unclear how this still-evolving, subsurface-drained basin would react. In regard to similar cases in the region, it can be expected that the lake will fill up, with subsequent rapid drainage.

4 Discussion

4.1 Glacier thermal regime and impact of glacier contact on lake development

The GPR data show a two-layered structure of the glacier, with prevailing cold ice and a temperate layer along the glacier bed. The terminus is frozen to the ground, and cold ice dominates the glacier base 450 m up-valley. At mid-length, on a narrow glacier snout around 3700 m a.s.l., the ice thickness increases to more than 50 m, and temperate ice appears at the glacier base. Temperate ice spreads along the base over the next 1 km and terminates below cold ice in the accumulation area. The cold–temperate transition along most of this section is located 30 to 50 m below the glacier's surface. The upper part of the glacier (above 3900 m a.s.l.) consists mainly of ice below the pressure-melting point and is thus connected to the permafrost environment in bedrock.

The presence of both types of ice confirms the polythermal regime of the Adygine glacier that was originally assumed only for large valley glaciers in the Tien Shan region (Dyurgerov et al., 1995). The distribution of temperate ice is similar to that of the characteristic polythermal glaciers described by Copland and Sharp (2001) and Etzelmüller and Hagen (2005). Following the classification scheme by thermal structure proposed by Irvine-Fynn et al. (2011), the Adygine glacier corresponds to a thermal configuration with the presence of temperate ice near the bed in the ablation area. The smaller extent of warm ice in the accumulation zone implies a different temperature structure than that of the Tuyuksu glacier (the Zailiyskiy Alatau), which has a maximum thickness of warm ice in the accumulation zone, where the cold–temperate transition reaches the glacier surface (Nosenko et al., 2016).

Lake 2 is a typical example of a lake formed when a glacier retreats over a rock step. The lake's growth was coupled with glacier recession until the terminus reached a higher elevation than the lake outflow. After the loss of contact, the lake stabilised and its formation ended. Lake 3 had a similar genesis and is currently in the stage of growth. Thanks to intensive thermal exchange between lake water and ice, a glacier terminus recedes by melting, and sometimes even by calving, in favour of lake area and depth. As a lake loses contact with a glacier tongue, some development can continue due to the melting of residual ice in the lake bed or dam. The extent of changes is, however, considerably reduced.

The case of Lake 1, formed in an intramorainic depression, demonstrates the difference in development. The lake's development was not terminated or stabilised after loss of contact with the glacier. Its further evolution was linked to buried ice, its exposure and melting, and the degradation of permafrost. The depression formed in the 1960s and filled with water two decades later; yet even today, the lake is subject to changes in basin morphology and drainage. Such lakes do not grow to large dimensions like typical moraine-dammed lakes (e.g. Tsho Rolpa and Dig Tsho); however, they can hold sufficient volumes of water to cause an outburst flood. In many cases, such lakes are even non-stationary, filling up only when a subsurface channel gets blocked, often in times of increased inflow due to snow or ice melting (Erokhin et al., 2017). The 2012 outburst from the Teztor Lake, situated in a neighbouring valley at a similar altitude to Adygine, serves as an example. When filled, the lake has a volume of 70 000 m³ and is usually drained by a subsurface channel with a flow rate of a few cubic metres per second (Erokhin et al., 2017). Nevertheless, in 2012, the channel's capacity enlarged, the flood changed into a debris flow, with peak discharge estimated at 350 m³ s⁻¹ (at the junction of the Teztor and Adygine valleys, ∼3 km from the glacier), and 200 000 m³ of debris was deposited at the fan entering Ala Archa Valley (Erokhin et al., 2017). The flood continued further downstream, causing dismay among the capital's inhabitants and material damages of ∼ USD 100 000 (Zaginaev, 2013). It is the often variable subsurface drainage, the common presence of buried glacier remnants, and a steep valley full of loose sediment downstream that make these lakes potentially dangerous.

4.2 Model uncertainties and future conditions

Huss et al. (2014) described, in detail, the uncertainties of glacier and runoff modelling in general. In our case, the GERM was calibrated using the summer and winter mass balance data of Golubina glacier from 1981 to 1994. Whereas simulated winter mass balances for Golubina showed no consistent trend of under- or overestimation, the modelled summer mass balance slightly overestimated the real conditions from 1981 to 1984 and generally underestimated it from 1985 to 1994. The performance of the model for simulating discharge is described by the Nash–Sutcliffe efficiency (Nash and Sutcliffe, 1970), which is 0.64 based on monthly mean values of the Adygine catchment for 1960–1987. The biases and uncertainties in the downscaled data for precipitation over the complex mountain topography remain a weakness and need further improvements in the future. Also, the lack of a DEM from 1960, as well as the glacier area from this time, attenuated the accuracy of the calibration of the GERM simulations.

Both the area and volume of glaciers all over Tien Shan are expected to decrease throughout the 21st century (Sorg et al., 2012). Aizen et al. (2007b) also predicted significant glacier degradation linked to the ELA shifting to higher altitudes because rising air temperatures will not be balanced out by sufficient increases in precipitation totals. Sorg et al. (2014) used the GERM for glacier and runoff modelling in Chon Kemin Valley (Zailiyskiy and Kungey Alatau, Tien Shan), where the glaciers are expected to have vanished by 2080 under the more pessimistic scenarios (dry–warm, wet–warm). The future runoff results of Chon Kemin (Sorg et al., 2014) are very similar to those of our study area, with “warm” scenarios expecting peak runoff in the 2020s and the main change in spring runoff (increase) being caused by higher winter precipitation and enhanced snowmelt in spring. A general prediction for runoff change in Tien Shan is a possible increase in the near future, followed by a steady decline until the end of the 21st century (Sorg et al., 2012). However, different regions in Tien Shan may exhibit varying responses and uncertainties connected to future runoff modelling. An advanced state of glacier degradation was presented by Huss and Fischer (2016), focusing on small-sized glaciers (<0.5 km²) in the Swiss Alps. Run-off from these glaciers has declined since the peak runoff year, which occurred between 1997 and 2004. There has been a significant decrease in August runoff, which is typical for the advanced stage of glacier runoff decline. Based on a comparison with expected future summer runoff from the Adygine glacier, we can estimate that Adygine is currently close to the stage of peak runoff.

Our approach to the identification of potential locations for future lake formation is based on glacier bed topography derived from in situ data. For obvious reasons, it is not possible to apply such a method to the whole region. Numerical models are a good substitute and are widely used in large areas, from the Peruvian Andes (Colonia et al., 2017) to the Swiss Alps (Linsbauer et al., 2012) and the Himalayas–Karakoram (Linsbauer et al., 2016). The commonly used model is GlabTop, introduced by Linsbauer et al. (2009) and Paul and Linsbauer (2012), with an automated version (GlabTop2) introduced by Frey et al. (2014). The model calculates glacier ice thickness, based on a DEM, and glacier outlines and branch lines; the automated version (GlabTop2) avoids manual delineation of the lines. A similar concept and results (Frey et al., 2014) were obtained from a model presented in Huss and Farinotti (2012). Potential sites with overdeepenings can also be identified, or confirmed, based on the three morphological criteria (glacier surface features) introduced by Frey et al. (2010). This is a more laborious approach, which includes manual analysis of the DEMs and high-resolution satellite/aerial imagery. If we apply these morphological criteria to the Adygine glacier, one location with a distinct narrowing, and also a change in slope, can be identified. This coincides with our results pointing to the presence of overdeepenings at the location of the current glacier tongue. Nevertheless, as Haeberli et al. (2016b) correctly pointed out, there are still significant limitations to understanding the principles of the depth erosion by glaciers and thus to modelling glacier bed overdeepenings. Moreover, the estimation of debris cover volume left after glacier recession is a challenge, as excessive debris can smooth topographic irregularities and lead to the formation of outwash plains, and not lakes (Linsbauer et al., 2016). We are also aware of the limitations to our approach. Due to the low resolution of the GPR-derived ice thickness raster, the resulting extent of potential sites for lake formation must be interpreted with caution.

4.3 Lake outburst susceptibility assessment

The presented evaluation strategy combined knowledge from the field with the benefits of remotely obtained data, providing a good insight into the interaction of the site's elements (a moraine, buried ice, perennially frozen ground, meltwater channels, a glacier, lakes) and a broader spatial and temporal perspective, thanks to the airborne imagery. Field investigations are necessary to produce a realistic assessment of lake outburst susceptibility. McKillop and Clague (2007) mentioned that knowledge of seepage and lake bathymetry and its changes is important information for understanding dam hydraulic conditions. Bolch et al. (2008) highlighted subsurface glacier meltwater routing, which we also believe to be extremely important, together with the buried ice distribution. At the Adygine ice-debris complex, the ERT technique helped us to discover an extensive ice core near Lake 2, buried under a considerably thick (8–10 m) layer of mostly frozen debris. A similar structure was observed, for example, at the moraine dam of the Thulagi glacier lake, Himalayas (Pant and Reynolds, 2000). The seepage routes and hydrological processes within moraine dams were also examined with SP and ERT in Thompson et al. (2012) at Miage glacier, Italian Alps. The resistivity values attributed to materials detected at our study site are consistent with the common resistivity value ranges of the different materials summarised in Kneisel (2006). Frozen ground can exhibit a wide range of resistivity values, varying from 10³ to 10⁶ Ω⋅m, depending on ice content, temperature (related to amount of unfrozen water), and content of impurities (Haeberli and Vonder Mühll, 1996). Glacier ice typically has relatively high resistivity values of well over 10⁶ Ω⋅m, which means it is practically non-conductive (Kneisel, 2003).

Remote sensing has proved to be a very useful tool in the field of environmental studies and is the basis for many hazard assessments (e.g. Allen et al., 2016; Bolch et al., 2008; Huggel et al., 2004; McKillop and Clague, 2007). Comprehensive hazard assessments elaborated as a basis for wide usage include, for example, the GAPHAZ Technical Guidance Document (GAPHAZ, 2017), the strategy prepared by ICIMOD, focused on the high-mountain region of the Himalayas–Karakoram (ICIMOD, 2011), and a useful summary of various lake outburst hazard assessments is also provided in Emmer and Vilimek (2013).

For a first, large-scale outburst hazard assessment, the trigger is often represented solely by a fall of mass into a lake (e.g. Allen et al., 2016). The presented procedure aims to encompass all trigger factors that could have the capacity to lower a dam's or lake basin's stability and cause lake drainage. However, there are still some limitations to a full understanding of the varying mechanisms of lake outbursts. We did not include permafrost degradation as a separate trigger, but we did consider it to be one of the main factors that can increase the outburst susceptibility in future, as presented for the Alps by Haeberli et al. (2017). The presence of permafrost at the site is supported by other studies (e.g. Marchenko et al., 2007), the global permafrost model of Gruber (2012), and the site's negative MAAT of recent years. Although the measurement period is rather short (2008–2013), the long-term MAAT of the site may, in fact, be even lower, given the low value (MAAT of −7.4 ^∘C) from another Kyrgyz station located at similar altitude (3614 m a.s.l.) in central Tien Shan (Bolch et al., 2019). An earthquake is definitely a phenomenon that can destabilise a structure via unconsolidated material, although McKillop and Clague (2007), for example, decided not to include it in their hazard assessment due to the low spatial variability. Several major earthquakes with epicentres in the vicinity of Kyrgyz Ala-Too were reported in the 20th century. The last case was the M7.5 quake near Suusamyr (∼50 km from Adygine) in 1992. Also, increased inflow due to heavy rainfall has been reported to cause moraine dam failure (Yamada, 1998), but such events are not common in the studied region. Heavy rains have been reported to trigger a debris flow (Zaginaev et al., 2016) but have not been linked to a lake outburst.

Glacier retreat in our site's case may not lead to new ice- or snow-avalanche-starting zones, but the rising MAAT could influence the glacier thermal regime and cause a shift in certain parts from cold- to warm-based, thereby lowering the friction between the glacier bed and steep ice masses (Huggel et al., 2008). Areas exposed after glacier recession often consist of unconsolidated material, which, if steep, may become a potential starting zone for landslides or debris flows (Haeberli et al., 2017). Already exposed steep slopes may become destabilised and fail due to the melting of interstitial ice or, where present, the degradation of mountain permafrost (Huggel et al., 2010). The formation of new lakes at higher elevations could increase the hazard of the site due to a possible cascade effect, where even a small volume release could trigger the outburst of a lower-lying lake. Such situations are expected to occur in the future, and the potentially arising challenges should be addressed in advance (Haeberli et al., 2016a).

5 Conclusions

This paper aimed to describe the impact of glacier retreat on its forefield, to highlight the differences in evolution and outburst susceptibility of proglacial lakes, and to promote a more field-based, regionally specified approach to hazard and susceptibility assessment.

The main results concerning the studied site are as follows:

• The glacier is subject to relatively rapid recession, which is comparable to that of other glaciers in Tien Shan. Since 1962, it has retreated, on average, by 11.5 m a⁻¹. By 2050, it is expected to have shrunk to 55.6 %–73.2 % of its 2012 extent.

• As a consequence, glacier runoff will change in the coming decades. Specifically, the beginning of the ablation season is likely to occur earlier, compared to the current situation, and the peak runoff from melting snow will very likely shift from mid-June to mid-May, resulting in a slightly earlier start to glacier ice melting.

• A three-level cascade of glacial lakes has formed in the glacier forefield, with the largest (Lake 2) retaining a volume of about 200 000 m³. However, due to its stable surface drainage and lack of contact with the glacier tongue, its development terminated in the 1990s. Two other lakes have the potential to grow – Lake 1, thanks to filling and the melting of buried ice, and Lake 3 by further glacier retreat. The latter is currently the most dynamic lake; having formed in 2005, to date, it has expanded to a volume of 106 000 m³ and a depth of >14 m.

• Lake 1 and Lake 3 were categorised as having medium outburst susceptibility, mainly due to their subsurface drainage and the presence of buried ice in their vicinities. Despite its large size, the outburst susceptibility of Lake 2 was assessed as low. The destabilisation of steep slopes, exposure and melting of buried ice, and changes in glacier runoff may further deteriorate the stability of Lake 1 and Lake 3, leading to an increase in outburst susceptibility in the future.

A very similar development in other glacier complexes in this region is expected, with runoff peaks occurring at present or in the near future, accompanied by permafrost degradation, the gradual melting of buried ice, and the formation of new lakes. This will result in new outburst hazards that could arise in a relatively short time. One of the main outputs of the study is an analysis of lake development dynamics, which indicate a lake's potential to grow and its related outburst susceptibility. As we have documented, lake size is not as important a factor, in terms of susceptibility, and we would rather emphasise the influence of buried ice in relation to the presence of perennially frozen ground and the unpredictability of subsurface drainage channels. In terms of future susceptibility change, we would like to draw attention to the influence of glacier runoff. The expected extension of the ablation season, and the related activity in glacial lakes, may lead to increased outburst susceptibility. The prolonged period of the melting season in relation to buried ice will result in a greater potential for the alteration of the hydrological network of drainage channels. Further research should be dedicated to the hydrological connections in proglacial morainic environments and the interactions of subsurface water flow with permafrost and buried ice, as there is still an insufficient understanding of these processes.