Between 2009 and 2020, 74 bathymetric surveys were conducted on 57 lakes (Figure 1; Supplementary Table S1). Of these, 13 lakes are managed and bathymetries of 10 lakes were measured repeatedly following lake level lowering. All surveys were conducted in late July—August when lake levels are at their maximum. Using a classification of lakes given in Kapitsa et al. (2017), we identified most of the surveyed lakes either as glacier-connected lakes (Type 1; 21 lakes) or glacier-disconnected lakes positioned on the moraines formed in the 20^th—21^st centuries (Type 2; 28 lakes including Lake Kapkan). Lake Kapkan (Jetysu Alatau) changed from the Type 1 lake in 2011 to Type 2 lake in 2018 and was included in both categories. Eight lakes were positioned on the Little Ice Age (LIA) or older moraines (Type 3), and one was a rock-dammed lake (Type 4). There were no ice-dammed and supra-glacial lakes in the study region.

Lake depths were measured using Lowrancе LMS-480 echo sounder with a dual 200 kHz and 50 kHz frequency transducer and a GPS receiver enabling simultaneous soundings and position measurements. The 200 kHz frequency was used throughout all measurements because of the relatively small depths of the lakes. Heave was not accounted for because the small size of the lakes and calm weather on the days of surveys ensured that the surface of the water was still. The vertical uncertainty of Lowrancе LMS-480 200 kHz channel was ±0.1 m. The GPS receiver had a high estimated horizontal uncertainty of ±9 m and to improve the accuracy of the measurements, the simultaneous position measurements were taken using Garmin GPSMAP 64 with an estimated accuracy of ±3 m.

The echo sounder was rigidly attached to a float fasten to an inflatable boat (Figure 2A). Lake depths were measured over the entire area of the lakes along the bathymetry survey lines spaced at approximately 10–20 m with the along-profile increment of 3–10 m depending on the size of a lake. If the signal indicated sharp changes in depth, the density of measurements was increased. The mean density of measurements varied between three and 30 soundings per 1,000 m² (Supplementary Table S1). In eight surveys (six lakes), the survey lines were spaced at about 50 m and these surveys can be identified in Supplementary Table S1 through the absence of the lake slope data (Figure 3 shows lake slope parameters).

The obtained data were imported to ArcGIS. Two methods of interpolation—triangular irregular network (TIN) and Kriging—were applied to the point data, and compared. The difference between the outcomes was within the uncertainty of measurements estimated using standard error propagation formula applied to two terms: i) vertical uncertainty of depth measurements (±0.1 m) times lake area and Garmin GPSMAP 64 accuracy (±3 m). TIN was selected as a more computationally efficient method. Raster data sets were generated from TIN; and bathymetric maps with depth isolines were generated from the raster data sets. Morphometric parameters of the lakes were calculated using ArcGIS functions and used in the compilation of the bathymetric inventory.

The inventory of the surveyed lakes consists of i) geographical coordinates and altitude (defined as altitude of water surface at a time of the survey); ii) type of lake (including glacier-connected with glacier tongue); iii) morphometric parameters derived from the bathymetric surveys (Supplementary Table S1); iv) type of dam and drainage (partially completed and not included in this paper); iv) details of lake management, past GLOF, and presence of lakes upstream and downstream (lake cascades). The obtained morphometric parameters are illustrated by Figure 3. Maximum length (referred to as length) is defined as the longest straight line between two points on the waterline. It is not related to the orientation of a valley or position of an outlet because not all lakes have visible outlets. Maximum width (referred to as width) is the longest distance along the axis perpendicular to the length of the lake. Lake elongation was calculated as a ratio between lake’s width and length. Mean depth was derived from the three-dimensional raster data sets and volume values were calculated bathymetrically using ArcGIS. Following Muñoz et al. (2020) and in order to facilitate the compilation of the consistent regional data sets, we calculated maximum and average front slope (MaxFS, AFS) and maximum and average back slope (MaxBS, ABS) of the lakes (Figure 3). Back slope refers to the side next to (or in the direction of) glacier tongue and front slope refers to the opposite direction. In both cases, the average slope was measured from the top of the water level to the maximum depth (Figure 3).

In addition to the data derived from the surveys, lake areas were measured using satellite imagery from July-September 1999, 2001, 2007, 2008, 2017 and 2021 (Supplementary Table S2) using methodology described in Kapitsa et al. (2017), Kapitsa et al. (2018).

To avoid bias introduced by the repeated surveys, duplicate measurements were not used in the statistical analyses (except Figure 4 and Supplementary Figure S1 later in the text). Data from the latest surveys were used for those lakes which were surveyed multiple times (Supplementary Table S1). Values of mean lake depths and volumes were regressed against other morphometric parameters of the lakes. Two goodness-of-fit metrics were used to assess regression models: coefficient of determination (R ²) and Mean Square Error (MSE). Following Huggel et al. (2004), relative error (RE) (defined as calculated minus measured value divided by calculated value) was used to compare measured and modelled depth and volume values. Standardized residuals (difference between measured and calculated values divided by standard deviation of the calculated value sample) were calculated and a standard threshold of ±2 was used to identify outliers. Statistical analyses were conducted using Minitab and R software.

Future hazard anticipation is becoming an important part of risk reduction. This approach relies on the prediction of future evolution of proglacial lakes. Here, we modelled the evolution of the existing glacier-connected Type 1 lakes using GlabTOP2 model which derives glacier bed topography based on glacier surface slope. GlabTOP2 is fully described in Frey et al. (2014) and its application to modelling the potential future lakes in the northern Tien Shan is discussed in Kapitsa et al. (2017). The model calculates ice thickness according to Eq. 1: h = τ f ∙ ρ ∙ g ∙ sin ⁡ α (1) where h is ice thickness m), τ is basal shear stress, f is shape factor (0.8), ρ is ice density (900 kg m^-3), g is acceleration due to gravity (9.81 m s^-2) and α is glacier surface slope. The basal shear stress is estimated from an empirical relationship between τ and glacier height extent Δz (i.e., maximum elevation minus elevation of the glacier tongue) according to Haeberli and Hoelzle (1995).

We used glacier outlines as in 2000 obtained from regular cataloguing of glaciers of Kazakhstan (Severskiy et al., 2016) and the Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM) surface topography to model overdeepenings in which glacier lakes may form and expand. To reduce the influence of small-scale surface undulations in the ice thickness modelling, SRTM DEM with 30 m resolution was down-sampled to 45 m cell size which has been shown to generate the optimal results in GlabTOP2 modelling (Kapitsa et al., 2017). The overdeepenings in glacier beds were identified following the methodology outlined in Linsbauer et al. (2012); Linsbauer et al., 2016 and Kapitsa et al. (2017).

Glaciers have retreated since the baseline year of 2000 in the study region and Type 1 glacier-connected lakes have formed and expanded (Kapitsa et al., 2017; Kapitsa et al., 2018). The morphometric characteristics of the overdeepenings, corresponding to the Type 1 glacier-connected lakes formed in the predicted locations including area, volume, maximum and mean depth were derived and compared to the measured values. A comparison between the measured and modelled parameters is indicative because most lakes continue growing and a discrepancy between the modelled and the observed values can be due to both, the model error and uncertainty associated with the further growth of the lakes. However, the repeated surveys showed that mean lake depth is a conservative parameter which changed with time less than area (which mostly controls increase in volume) and we focused on assessing the ability of GlabTOP2 to simulate mean depth and bathymetry of the overdeepenings.

Table 1 and Figure 4 present the descriptive statistics of morphometric parameters of the surveyed lakes. The lake areas are relatively small ranging between 23.39x10⁴ (Maksimov Lake in the Kungey Alatau; Figure 2B; Figure 5A) to 0.2x10⁴ m² but their median (2.87x10⁴ m²) and mean (5.55x10⁴ m²) areas are comparable with a mean area of 2.7x10⁴ m² derived from a larger sample of over 600 lakes in the Jetysu Alatau (Kapitsa et al., 2017). The Type 1 (glacier-connected) lakes tend to have larger areas than the Type 2 lakes (forming on the young moraines) with the median areas of 5.75x10⁴ and 2.0x10⁴ m², and mean areas of 6.85x10⁴ and 4.0x10⁴ m², respectively. The Type 3 lakes (forming on the LIA moraines) had a mean area of 6.28x10⁴ m². The sample include eight lakes only and their larger area is an artefact of focusing on larger lakes at risk of GLOF which are otherwise infrequent among the Type 3 lakes.

The surveyed lakes are relatively shallow. The median values of mean and maximum depths are higher for the Type 1 lakes at 8.3 and 19.0 m, respectively, against 4.6 and 10.5 m for the Type 2 and 3 lakes. The deepest lake in the sample is the glacier-connected Maksimov Lake which has the mean and maximum depths of 22.9 and 57.6 m, respectively (Figure 2B; Figure 5A). The deepest parts of the Type 1 glacier-connected lakes are usually those adjacent to the glacier tongues while in the Type 2 lakes, the middle parts of the lakes tend to be deeper. The average and maximum slope values for the surveyed sample are shown in Table 1. The median and mean values of ABS of Type 1 lakes are 18^o and 29^o against the respective AFS values of 13^o and 18^o while the median MaxBS value is 70^o. For Type 2 lakes, the median and mean ABS values are 13^o and 14^o and the median value of the MaxBS is 32^o. There is no significant difference between the AFS and MaxFS values of Type 1 and Type 2 lakes. Overall, the surveyed lakes have simple bathymetries due to their small size.

The largest median values of volume characterise the Type 3 lakes, which have larger areas and smaller depth, followed by Type 1 lakes which are deeper (Figure 4). The Maksimov Lake (Type 1 lake which is currently not managed) has the largest volume of 53.89x10⁵ m³.

Values of mean lake depths and volumes were regressed against other morphometric parameters of the lakes using data from the latest surveys. Regression against lake area, length, and width produced coefficients of determination significant at 5% confidence level (Supplementary Figure S1; Table 2). The relationships were stronger for the Type 1 lakes for which linear and polynomial relationships produced better results than the frequently used logarithmic and power fits with the R ² values of 0.76–0.78 and MSE of 7 m for the depth-area relationships (Table 2). For the Type 2 lakes, both lake area and width exhibited the strongest correlation with depth with R ² of 0.50–0.63 and MSE of 5–7 m (Table 2). The low MSE values are low implying small bias in the simulation of lake depths. For the whole sample, lake width and area were better predictors of depth than lake length but R ² values were lower than for the Type 1 group (Table 2) which has more homogenous characteristics (Figure 4). Other parameters did not correlate with depth either for the whole sample or for the Type 1 and Type 2 sub-samples. In particular, lake elongation, which has been used in other studies as a predictor of depth and volume (Qi et al., 2022), did not show significant correlation with either depth or volume in this data set.

Two lakes, Boskol and Lake N 4 in the Aksu catchment in the Jetysu Alatau (Supplementary Table S1) were identified as outliers in the regressions between mean depth and other morphometric parameters. Lake Boskol is a Type 3 lake formed on the LIA moraine. It has the third largest area of 19.75x10⁴ m² but its depth does not exceed 7 m. Lake N 4 has usual shape with a length of 1,050 m (the largest in the data set), width of 230 m and elongation of 0.22.

While relationships between lake volumes and other morphometric parameters, especially area, are criticized for the inherent autocorrelation, they are useful in practical applications (Huggel et al., 2002; Wang et al., 2012; Loriaux and Casassa, 2013; Cook and Quincey, 2015; Khanal et al., 2015; Wood et al., 2021; Medeu et al., 2022). Results of regression between the bathymetrically derived volumes and other morphometric parameters are shown in Table 3. The highest coefficients of determination were obtained for the Type 1 glacier-connected lakes with R ² values ranging between 0.71 and 0.99. Lake area and mean depth were the best predictors. Regressing volume against mean depth achieved simultaneously the highest R ² and the lowest MSE values for the glacier-connected Type 1 lakes and for the whole sample. In general, power equations show the best fit but polynomial functions can be used to strengthen these further, e.g., between volume and area of Type 1 lakes where the use of polynomial function produces lower MSE value. There was substantial uncertainty in the estimation of volumes of larger lakes from depth and area when using power equations. Thus, volumes of the three largest lakes—Maksimov, Kapkan and Taskol with areas around 20x10⁴ m² and volumes in excess of 30x10⁵ m³—were underestimated by 30%–40% (Figure 6A). The second-order polynomial equation provides a more accurate approximation of the volumes of Type 1 lakes including those of the larger lakes (Table 3; Figure 6A).

The ability of the morphometric parameters to predict volume was lower for the Type 2 lakes and the combined data set containing all lakes. Using area as predictor allowed maximizing R ² and minimizing MSE for the Type 2 lakes. In addition to area, maximum depth and lake width were the two parameters which predicted volume of Type 2 lakes best. There was no difference between the performance of the polynomial and power equations (Figure 6B).

The obtained power equations linking mean depth and area (Table 2) and volume and area (Table 3) were tested using the independent data sets of lake depths and bathymetrically derived volumes compiled by Huggel et al. (2002), Wang et al. (2013), Cook and Quincey (2015), Zhang et al. (2023) and the original data by Janský et al. (2010), Muñoz et al. (2020), and Zhang et al. (2023). Duplicate measurements were excluded from all data sets. Data for volumes only were used from Zhang et al. (2023) because maximum instead of mean depth was reported in this study. Two lakes [Leones in Patagonia (Harrison et al. (2008); Loriaux and Casassa (2013)] and Galong Co [China (Zhang et al., 2023)] were excluded from the assessment of volume calculations because their volumes are an order of magnitude higher than those of any other lake in the sample (Supplementary Table S3). Data for the lakes with natural dams only were used from Muñoz et al. (2020). The compiled data includes three types of lakes (Table 4), two of which—ice-dammed and supra-glacial—did not occur in our study area. We could not verify the accuracy of the published measurements excluding only the data for the Brazo Rico Lake in South Patagonia from the Cook and Quincey (2015) data set because depth measurements were provided for the lake areas adjacent to the glacier tongue instead of the mean depth (Stuefer et al. 1999; Stuefer et al., 2007) and were a clear outlier. The combined samples contained 159 and 188 lakes for the assesment of depth and volume calculations, respectively (Supplementary Table S3). On average, lakes in the published data sets were larger than in our bathymetrical data set with means and standard deviations of (5.91 ± 10.26)x10⁵ against (5.6 ± 5.9)x10⁴ m² for areas and 21.6 ± 23.1 and 6.8 ± 4.7 m for depth, respectively.

Table 4 shows the results of the application of the power area-depth and area-volume equations for the Type 1 glacier-connected lakes (Tables 2 and 3) which showed marginally better performance than the equations for the Type 2 glcier disconnected and all surveyed lakes. There was a close correlation between the modelled and measured depths of the moraine, rock and ice-dammed lakes with the R ² value of 0.62—0.68 with and R ² value of 0.64 for all lakes (Figure 7A). On average, the modelled values were 5–7 m lower than the measured values. The errors of the regression models, however, were higher with the lowest MSE value of 190 m obtained for the moraine and rock-dammed lakes. A high MSE value of 670 m obtained for the ice-dammed lakes results from a large error in the calculation of mean depth of Lake No Lake (155 m against the calculated value of 56 m). Removal of this lake from the assessment reduces MSE for the ice-dammed lakes to 80 m. By contrast, depth of the supraglacial lakes was overestimated. Correlation between the measured and modelled depths of the supraglacial lakes was strong (R ² of 0.87) but it was derived from a very small sample.

RE were calculated to assess the ability of the new equations to predict depths and volumes of individual lakes. The lakes were binned into four categories according to the assigned RE (Table 5). Cook and Quincey (2005) characterised lakes whose depth was calculated with an error of over ±49% as “unpredictable”. A large discrepancy, however, can reflect systematic bias. In this study, standardized residuals exceeding a threshold of ±2 were used to identify outliers.

Overall, depths of 55% of all lakes were calculated with an error less than ±49% (Table 5). The standardized residuals were generally higher for lakes with depth greater than 20 m but the critical threshold was exceeded for seven lakes only. The highest underestimation characterised two large lakes—Nef in North Patagonia (Loriaux and Casassa, 2013) and Lake No Lake in British Columbia (Geertsema and Clague, 2005)—with areas greater than 50x10⁵ m² and depth of 150 m. Both lakes are deep for their surface areas. Other formulae were applied to and failed to predict the depth of Lake No Lake (Cook and Quincey, 2015). On average, the lowest error (±34%) was obtained for the ice-dammed lakes exceeding 49% only for Lake No Lake. The sample of ice-dammed lakes is relatively small (Table 4) which makes its sensitive to outliers. If Lake No Lake is removed from the sample, R ² between the measured and modelled depths of the ice-dammed lakes increases to 0.84 and RMSE changes from ±25.9 m to ±9.3 m. Predictability of the moraine- and rock-dammed lake depth varied strongly within the sample. In addition to Lake Nef, depths of lakes Maud in New Zealand (Allen et al., 2009) and Auquiscocha, Purhuay, and Rusgo in the Andes (Muñoz et al., 2020), exceeding 40 m, were strongly underestimated. However, correlation between the measured and modelled depths was stable and not significantly affected by the removal of the outliers. Among the supraglacial lakes, a shallow Lake Petrov in the Tien Shan is an outlier whose depth was not reliably predicted by the previously published scaling equations (Cook and Quincey, 2015). The lake is positioned on a debris-covered glacier terminus and experiences strong seasonal fluctuations in area and depth (Janský et al., 2010).

There was a strong correlation between the calculated and bathymetrically measured volumes of lakes with R ² value of 0.80 for the overall sample (Table 4; Figure 7B). Volumes were calculated with an error of less than ±49% for 44% of the lakes (Table 5). Derived from a sample of smaller lakes, the scaling equation generally underestimated volumes except for the supraglacial lakes whose volumes were overestimated (Table 4). Results for the ice-dammed lakes were affected by the inclusion Lake No Lake in the sample and improve significantly if this outlier is removed with R ² value increasing to 0.94 and a smaller difference between measured and calculated values which became 96.05x10⁶ and 73.53x10⁶ m³, respectively. In addition to Lake No Lake, four other large lakes—Nef (770.7x10⁶ m³) in Patagonia (Loriaux and Casassa, 2013), Tasman (510.0x10⁶ m³) in New Zealand (Robertson et al., 2012), Phantom (500.0x10⁶ m³) in Axel Heiberg Island (Maag, 1963), and Rewoco (139.0x10⁶ m³) in China (Zhang et al., 2023)—were identified as outliers (with standardized residuals exceeding the threshold of ±2) whose volumes were significantly underestimated (Figure 7B). The Petrov Lake in the Tien Shan Janský et al. (2009) was an outlier whose volume was overestimated by 61%.

To expand the samples of measured depths and bathymetrically derived volumes of large lakes and those types of lakes which do not occur in the study area, we combined the newly obtained and published data. The measured depth and bathymetrically derived volume data sets included data for 216 and 245 lakes with natural dams excluding duplicate measurements (Figure 8A; Supplementary Table S3). There is a good correlation between lake depth and area with R ² values of 0.67 overall and 0.71 for the moraine- and/or rock-dammed lakes. This is comparable with the results obtained in this study for the more homogeneous group of smaller glacier-connected lakes (Table 2) and a significant improvement on the results by Cook and Quincey (2015) who reported stronger variability in lake depth for any given area with R ² of 0.38 derived from a smaller data set of 42 lakes. Analysis of standardized residuals shows that three lakes remain outliers in the combined depth data set (Figure 8A) and their removal from the data set reduces MSE by half. The values of standardized residuals are higher for the lakes with depths of 20–60 m than for lakes with depths of under 20 m.

The obtained bathymetries provided a unique opportunity to evaluate the skill of GlabTOP2 in simulating bathymetries of 18 Type 1 glacier-connected lakes. The measured lake areas and volumes in this sample varied between 0.79x10⁴ and 23.39x10⁴ m² and between 0.22x10⁵ and 53.89x10⁵ m³, respectively. The modelled and measured mean depth values (±standard deviation) were 13.1 ± 5.9 and 9.3 ± 5.6 m, respectively. The mean maximum depth values were 32.5 ± 15.3 and 21.3 ± 12.5 m, respectively. Lake 13 in the Bien catchment is an outlier (discussed further in the text) and excluded from this comparison.

Figure 5 further illustrates the observed and predicted evolution of the selected lakes identified as dangerous (Kapitsa et al., 2018; Medeu et al., 2022). There is a good agreement between the modelled and measured bathymetries of the Maksimov Lake (Figure 5A). The measured and modelled maximum and mean depths values were 58 and 22 m against 50 and 20 m, respectively. The lake’s area and volume are projected to reach an of 45.45x10⁴ m² and 107.39x10⁵ m³, respectively, and a small slope of the glacier tongue (Figure 2B) indicates that the lake will indeed continue growing. A single bathymetric survey of the lake was conducted in 2018 and there is no information about actual changes in its depth. The multiple measurements of its area using satellite data showed that it doubled between 1999 and 2021 although the rate of area growth slowed down from 5% per year in 1999–2007 to 2.5% per year in 2017–2021 (Figure 5A). This comparison provides confidence in and raises concerns about the predicted growth of the lake which is projected to remain the largest in the Ile-Kungey Alatau. Although the lake’s dam is composed mostly of solid rock reducing the risk of GLOF, the lake is partially contained by a moraine dam (Figure 2B). Seismic activity, common in the region, may affect the rock dam too raising level of the potential hazard which is predetermined by lake volume and maximum peak discharge controlled by volume (Emmer and Vilímek, 2013).

There is a close agreement between the modelled and measured depth of Lake N 2 in the Turgen basin (Figure 5B). Its area increased from 2.25x10⁴ to 8.05x10⁴ m² between 1999 and 2021, however, GlabTOP2 does not predict any further increase. The Makarevich Lake in the Kaskelen basin is increasing in line with the GlabTOP2 results although the modelled depth in the sector of the lake adjacent to the glacier exceeds the measured values (Figure 5C). While the model error is not excluded, the regular drainage of the lake via the artificial channels may contribute to this discrepancy. Repeated surveys of the Galina Lake in the Uzun-Kargaly basin are required to establish whether the lake’s depth, currently reaching 19 m, is likely to reach the projected depth in excess of 30 m (Figure 5D).

Lake 13-bis in the Ulken Almaty catchment was identified one of the most dangerous lakes in the Ile Alatau (Kapitsa et al., 2018; Medeu et al., 2022). Its area and volumes increased from 3.05x10⁴ to 4.3x10⁴ m² and 1.8x10⁵ to 3.05x10⁵ m³, respectively, between 2011 and 2015 but declined afterwards due to the implementation of artificial lake drainage. Its mean and maximum depths reached 7 and 17 m, respectively, before 2015. These values are consistent with the GlabTOP2 predictions for the same sector of the lake (Figure 5E). Further increase in the lake’s area and depth are projected and further bathymetric surveys are required to establish if the lake is likely to reach the projected depths. By contrast, no further increase is projected for one of the most dangerous and the second largest in the data set Lake Kapkan (37.25x10⁵ m³) in the Jetysu Alatau (Kapitsa et al., 2017). Observations confirmed that the lake changed from Type 1 in 2011 to Type 2 in 2018 (Figure 5F) and there is no further potential for increase in lake area and volume.

Overall, GlabTOP2 simulates realistic depth values although there is a tendency towards the overestimation in comparison with the measured depths. In the assessed data set, there was a single outlier—Lake 13 in the Bien River catchment, Jetysu Alatau—where the parameters of the modelled overdeepening substantially exceeded the measured parameters (Figure 5G). The projected maximum and mean depth were 80 and 39 m, respectively. The maximum and mean lake depths measured in 2018, were 10 and 5 m in the sector of the lake where the largest depths were modelled. Importantly, the whole overdeepening, which now contains two lakes (Figure 5G), does not reach the depth and volume projected by the model.

Information about volumes of glacier lakes is important for modelling GLOF and runout distances of the floods. However, there are few direct measurements of lake depth and volumes which can only be obtained through bathymetric surveys. The current data set, derived from 74 surveys of 57 lakes, significantly improves data availability. It enables detailed characterization of lake bathymetries improving risk assessments by the regional risk-reduction agency KSMP in the Tien Shan where GLOF events cause substantial damage to the downstream communities and ecosystems (Medeu et al., 2022). It also enables better understanding of statistical relationships between morphometric parameters of glacial lakes in the region. It helps to establish more robust empirical relationships between area, depth and volume of glacial lakes globally when combined with the previously published data.

Area-depth and area-volume relationships, developed from the regional and global data sets, were widely applied in Central Asia (Konovalov, 2009; Konovalov and Rudakov, 2016; Medeu et al., 2022). Table 7 and Supplementary Figure S2 show a comparison between the bathymetrically derived values of lake volumes and depths and values calculated from the lake areas using equations published in the literature and derived in this study. Most equations produced similar R ² values. There was a closer correlation between the estimated and measured values of lake volumes (R ² values of 0.78–0.97) while there were larger discrepancies between the calculated and measured depths (R ² values of 0.50–0.68). For both parameters, there was a closer agreement with the measured parameters of the Type 1 glacier-connected lakes. The mean absolute values of RE were mostly moderate ranging between 25% and 40% and were lower for the Type 1 lake depth (Table 7). Higher errors in the calculation of volumes of the glacier-connected lakes resulted from the formula derived by Wang et al. (2012) based on the published surveys of the moraine-dammed lakes in the Himalayas. The application of the formula proposed by Fujita et al. (2013) for the calculation of lake depth in the Himalayas resulted in a strong overestimation of lake depths probably because it was derived from a sample of much larger lakes.

The area-depth relationships by Huggel et al. (2002), Wang et al. (2012) and the power equation for Type 1 lakes derived in this study produced the lowest mean absolute RE of 26% for lake depth. The formula by Khanal et al. (2015) produced the lowest mean absolute RE in the estimation of lake volumes with the formulae by Cook and Quincy (2015), Huggel (2002) and the one derived in this study producing similar results. Analysis of the actual RE values showed that most formulae overestimated depth and volume of smaller lakes with mean depths less than approximately 15 m (Supplementary Figure S2) generating positive bias in the calculated values (Table 7). The newly developed scaling equations do not have this limitation as errors of the depth and volume regression models are randomly distributed while the use of polynomial functions helps to avoid underestimation of depths and volumes of larger lakes. We note, however, that equation developed by Evans (1986) also produced randomly distributed errors and no systematic bias in the calculated depths and volumes of lakes with depths under 15 m (Supplementary Figure S2) while the equation by Popov (1991) produced randomly distributed errors for volumes of the shallow lakes.

All power functions underestimated depths and volumes of the lakes with mean depths and volumes exceeding 15 m and 3x10⁶ m³, respectively. (Supplementary Figure S2). The formula by Popov (1991) produced the largest errors underestimating volumes of larger lakes by 38%–82%. The use of the second-order polynomial function reduced the error of simulating volumes of larger lakes to 0%–4% and is recommended for a more accurate approximation of depths and volumes of the larger glacier-contact lakes in the study region.

In the global context, the main contribution is in a significant increase of the available bathymetries of glacial lakes. This helps constrain uncertainty of the application of area-depth and area-volume equations to the large and non-homogenous data sets and enables researchers to devlop sub-samples of lake depths and volumes to derive scaling equations for lakes of different types and sizes. The power equations derived from the combined data set of over 260 lakes (Supplementary Table S3) had the R ² values of 0.67–0.71 and 0.80–0.85 for depth and volume, respectively (Table 6). This is a considerable improvement on the results derived by Cook and Quincey (2015) from a smaller global data set. However, uncertainty in simulating lake depths and volumes using the combined data set increases with lakes’ size. While the values of standardized residuals do not exceed ±1 for the lakes with depth below approximately 20 m, there is a notable trend for higher residuals for lakes deeper than 20—30 m. Therefore, more measurements of larger lakes, which tend to be deeper, are needed. Supraglacial lakes form a distinct sub-sample of lakes whose depths are consistently underestimated by the developed equations (Table 4). This limitation was also noted by Cook and Quincey (2015). There is a limited number of surveys of supraglacial lakes globally (Table 4, Supplementary Table S3) and wider data acquisition is required.

The results of the bathymetric and remote sensing surveys are used by practitioners (KSAMP) in conjunctions with the GlabTOP2 projections (Kapitsa et al., 2017; Kapitsa et al., 2018) to improve efficiency of lake monitoring and management by focusing on those lakes which are reliably predicted to grow in the future in addition to those posing immediate risk. Lake management involves lowering lake levels by installing syphons and construction of the artificial drainage channels. This practice proved efficient in the prevention of GLOF whose frequency peaked in the 1970s when strong positive temperature anomalies and glacier melt were observed in the region (Shahgedanova et al., 2018). In 1971–1980, 20 GLOF occurred while a single lake was managed. By contrast, there were four GLOF since 2001 despite the observed climatic warming and accelerating glacier retreat because KSMP significantly increased the number of managed lakes to 26 (Supplementary Figure S3) reducing risk of GLOF which was otherwise identified as high on the global scale (Taylor et al., 2023). However, lake level lowering measures are expensive and require the deployment of extensive workforce and machinery in the challenging mountain environments (Medeu et al., 2022). Bathymetric surveys and modelling of lake expansion are important tools enabling decision-makers to optimize the risk reduction measures.

The repeated bathymetric surveys of Lake Kapkan (37.25x10⁵ m³) in the Jetysu Alatau, which is the second largest in the data set and one of the most dangerous in the region (Kapitsa et al., 2017), showed that its mean and maximum depths were reduced from 19 to 41 m, respectively, in 2011 to 14 and 36 m in 2018 following the artificial lake drainage. A significant reduction in the volumes of the extremely dangerous lake 13-bis in the Ulken Almaty catchment and Lake Manshuk in the Kishi Almaty catchment have been achieved since 2015 (Supplementary Table S1). Monitoring of Maksimov Lake was initiated on the basis of the GlabTOP2 modelling output which predicted a strong growth of the lake (Kapitsa et al., 2018) while the projected growth of Makarevich Lake provided justification for funding of the lake level lowering measures. The KSAMP aims to implement risk avoidance measures in addition to risk reduction by advising regional authorities and planning agencies on future dangers of lake formation in specified locations. This approach has become a part of the regional strategy of adaptation to climate change. It is, therefore, important that the model outputs are validated.

The comparison of measured and modelled bathymetries is indicative because most lakes continue growing and it is not possible to use any quantitative metrics. The inspection of the measured and modelled bathymetries of the Type 1 glacier-connected lakes (Figure 5) showed that most simulations provided realistic results although it appears that GlabTOP2 tends to overestimate lake depth. In the case of Lake 13 in the Bien catchment (Figure 5G), ClabTOP2 failed to reproduce the measured parameters of the overdeepening containing the lake. We attribute it to the significant overestimation of ice thickness: the simulated maximum depth was 168 m against the typical measured values of 30–60 m (Cherkasov, 2004).

We attempted to assess uncertainty by comparing the modelled ice thickness with ice thickness values compiled by Farinotti et al. (2019) from four different models: Huss and Farinotti (2012), Frey et al. (2014) (who used GlabTOP2 but with settings different to the ones applied in this study), Maussion et al. (2019), and Fürst et al. (2017). The values of maximum ice depth and its standard deviation were derived for each glacier terminating in the lake for which overdeepenings were modelled and measured bathymetric data were available. Overdeepenings were derived from each of the four models’ ice thickness data sets and their bathymetries were compared with our simulations using GlabTOP2 and with measured bathymetries (Supplementary Table S4). Overall, 15 glaciers and lakes were used in this comparison because data in the Farinotti et al. (2019) data set were not available for all glaciers. In case of managed lakes, we used bathymetries obtained before lake level lowering commenced.

This comparison has several limitations. Firstly, in our study, glacier outlines from 2000 (or nearby year) were used as contemporaneous with the SRTM data acquisition. The published ice thickness values were based on the glacier outlines from the Randolph Glacier Inventory (RGI) 6.0 (RGI Consortium, 2017) obtained in later years. This difference affects modelled bathymetries of the overdeepenings under the rapid glacier retreat in the study region. For example, one of the most dangerous Lake 13-bis in the Ulken Almaty catchment began to form in 2000 while by 2015, its volume had reached 4.3x10⁵ m³ and lake lowering measures commenced (Supplementary Table S2). Therefore, our simulation generates a larger area of its overdeepening than the simulations based on the later glacier outlines. Secondly, the assessed sample of lakes was small while the between-model variability was high (Supplementary Table S4).

Despite the limitations, the comparison of different ice thickness simulations highlighted indicators of errors in modelling overdeepenings using GlabTOP2. The mean values and standard deviations of ice thickness of the selected 15 glaciers were 73 ± 15 m in our study and 96 ± 19, 75 ± 16, 84 ± 20, and 61 ± 14 for models one to four, respectively. A strong between-model discrepancy in the simulated ice thickness was observed in a number of cases and was indicative of error in modelling overdeepenings. In case of Lake Bien, the maximum ice thickness of the adjacent glacier varied between 164 m (Huss and Farinotti, 2012), which is close to the value obtained in this study, and 99 m (Fürst et al., 2017) and all models appeared to have overestimated ice thickness (Cherkasov, 2004). In case of Maksimov Lake, which was simulated accurately using the GlabTOP2 data (both, from this study and from model 2 (Frey et al., 2014)), there was a closer agreement between all models in simulation of the maximum ice thickness with values ranging between 100 m (Fürst et al., 2017) and 130 m (Frey et al., 2014). Closer examination of the simulated ice thickness showed that overestimation occurred where glaciers change directions of their flow. These two criteria can be used in the future in the assessment of uncertainty.

Of all assessed models, GlabTOP2, as used in this study and by Frey et al. (2014), provided realistic simulations of overdeepening bathymetries (Supplementary Table S4; Supplementary Figures S4–S6). Results based on ice thickness from Frey et al. (2014) underestimate lake depth predicting the mean depth of future overdeepenings as 7.0 m against the currently observed value of 8.4 m (Supplementary Figure S4). For example, the Maksimov Lake was predicted to achieve mean depth of 12.6 m against the currently observed 22.9 m (Supplementary Figure S5). Other three models significantly underestimated area, volume, and depth (with average values of future mean depth of 2.0—3.6 m) and are not recommended for the simulation of overdeepenings in the study area. We note, however, that model 1 (Huss and Farinotti, 2012) provided a more realistic representation of bathymetry of Lake Bien 13 (Figure 7G) which was overestimated by both GlabTOP2 simulations (Supplementary Figure S6).