Glacier Ice Thickness Estimation and Future Lake Formation in Swiss Southwestern Alps—The Upper Rhône Catchment: A VOLTA Application

Abstract

Glacial lake formations are currently being observed in the majority of glacierized mountains in the world. Given the ongoing climate change and population increase, studying glacier ice thickness and bed topography is a necessity for understanding the erosive power of glacier activity in the past and lake formation in the future. This study uses the available information to predict potential sites for future lake formation in the Upper Rhône catchment located in the Southwestern Swiss Alps. The study integrates the latest available glacier outlines and high-quality digital elevation models (DEMs) into the Volume and Topography Automation (VOLTA) model to estimate ice thickness within the extent of the glaciers. Unlike the previous ice thickness models, VOLTA calculates ice thickness distribution based on automatically-derived centerlines, while optimizing the model by including the valley side drag parameter in the force equation. In this study, a total ice volume of 37.17 ± 12.26 km³ (1σ) was estimated for the Upper Rhône catchment. The comparison of VOLTA performance indicates a stronger relationship between measured and predicted bedrock, confirming the less variability in VOLTA’s results (r² ≈ 0.92) than Glacier Bed Topography (GlabTop) (r² ≈ 0.82). Overall, the mean percentage of ice thickness error for all measured profiles in the Upper Rhône catchment is around ±22%, of which 28 out of 42 glaciers are underestimated. By incorporating the vertical accuracy of free-ice DEM, we could identify 171 overdeepenings. Among them, 100 sites have a high potential for future lake formation based on four morphological criteria. The visual evaluation of deglaciated areas also supports the robustness of the presented methodology, as 11 water bodies were already formed within the predicted overdeepenings. In the wake of changing global climate, such results highlight the importance of combined datasets and parameters for projecting the future glacial landscapes. The timely information on future glacial lake formation can equip planners with essential knowledge, not only for managing water resources and hazards, but also for understanding glacier dynamics, catchment ecology, and landscape evolution of high-mountain regions.

1. Introduction

High-mountain regions are dominated by glacial and periglacial processes which are mainly governed by climate conditions [1]. Currently, the global temperature is ~1.5 °C above the pre-industrial level and is continuing to rise with a predicted range of 0.8 °C to 1.2 °C between the 2030s and 2050s [2]. As a consequence, a potential ice volume loss of ~65–80% between the early 21st century and 2100 under moderate warming and a near-complete loss under warmer conditions are expected for the entire European Alps [3,4,5,6].

One of the widely observed signs of climate change and ice loss has been the transformation of glaciers into a new landscape of rocks, debris, ice-debris complexes, and most importantly glacial lakes [7,8,9,10,11]. The combination of these features can form a highly dynamic and hazardous condition for the downslope settlements (e.g., Glacier Lake Outburst Floods (GLOFs)) [12,13,14,15,16,17]. The formation of proglacial lakes can also affect the stability of mountain glaciers and can partly disengage glacier behavior from climatic perturbations, with such processes having occurred during previous phases of deglaciation [18]. Despite their hazardous nature, the newly emerging proglacial lakes and ice-free basins have the potential to serve as geo-tourism attraction sites or water reservoirs for energy production [19,20,21,22]. Recently, understanding lake formation and modeling the potential location of future lakes has become a necessity not only for timely water resource and hazard management in high mountain regions, but also for glaciological and ecological research concerning the reconstruction of ice motion, landscape evolution [18,23,24], and high-mountain biodiversity [25,26,27,28,29,30,31,32], respectively.

Enhancements in the spatial resolution of satellite images and the advent of Unmanned Aerial Vehicles (UAVs) [33,34,35] have greatly improved the workflow of geomorphometry-related approaches [12,36], the production of digital elevation models (DEMs), and glacier mapping [37,38,39,40,41,42]. All of these improvements have finally resulted in distributed ice thickness models [43,44] for remote mountain environments. The primary application of such detailed approaches in lake formation studies is the projection of future ice-free mountain landscapes, particularly where the ice-water-sediment feedback has largely contributed to the formation of subglacial basins [7,24,45].

Currently, a wealth of approaches are available in this context. Among them, the scalar approach is the simplest method that has been widely used in glaciological studies. Such approaches only estimate the mean ice thickness and total volume of a glacier assuming that larger glaciers are thicker. However, several same-area glaciers have previously shown significantly different ice thickness values within the same geographical setting. Apart from this, a single mean ice thickness value is of little use for predicting lake formation [46,47]. Later, the methods improved based on unified glacier mapping approaches, better resolution of satellite images, and the availability of digital terrain data. These datasets provided topographic information for the detection of glacier outlines, elevation ranges, and slope within the glacier boundaries. Through these advancements, methods that yield distributed ice thickness and ice volume emerged (e.g., slope-dependent approaches based on basal shear stress and perfect plastic behavior theories) [37,43,48,49,50]. Over the past decade, more complicated models have become increasingly available. Such approaches consider additional datasets such as mass conservation, ice flow dynamics, surface velocities, and mass balances [51,52,53,54]. Further, machine learning techniques for simulating glacier beds from deglaciated areas [55], and Bayesian inference [56], paved the path for the next generation of ice thickness modeling and automation.

This study employs the Volume and Topography Automation (VOLTA) model [57], a perfect-plasticity based model, to estimate distributed ice thickness for the glaciers of the Upper Rhône catchment. VOLTA has the advantage of considering the effect of side drag variations in glacier margins on ice thickness estimation [49]. In addition, it fully automates glacier flowline delineation and calculates basal shear stress for each individual glacier, thus making it more effective for regional studies. Previously, VOLTA has been successfully applied to southern South America [58], and the entire Antarctic Peninsula [59]. Further, a modified version (REVOLTA) was introduced for the reconstruction of glaciers at the Last Glacial Maximum across the New Zealand Southern Alps [60]. However, to date, the model has not been applied to the Swiss Alps for the purpose of future lake prediction. Against this background, we incorporate the latest available Swiss Glacier Inventory (SGI) and a high-resolution DEM in VOLTA and assess the capability of this model as an alternative to the previously used models. We also develop a simple framework for delineating potential lake formation sites based on VOLTA ice thickness results, and morphometric criteria suggested by Frey et al. [61] and refined by Magnin et al. [62]. Accordingly, the objectives of this study can be summarized as: (i) estimation of distributed ice thickness within the Upper Rhône catchment, (ii) modeling subglacial bed topography, (iii) prediction of future lake sites, (iv) investigation of geometry and morphological characteristics of predicted lakes, and v) validation of results with field measurements and comparison of the VOLTA performance against Ice Thickness Estimation Method (ITEM), Glacier Bed Topography (GlabTop) [50] and GlabTop2 [63] models.

2. Study Area

The Upper Rhône catchment is located in Canton Valais in southwestern Switzerland (Figure 1), covering a total area of ~5290 km² with 340,000 inhabitants [64]. Enclosed by two alpine mountain ranges, the Bernese Alps in the north and the Pennine/Valais Alps in the south, the catchment elevation range extends from 373 m to 4631 m, with a mean altitude ~2085 m. The north side of the Bernese Alps is notably wet under oceanic conditions, while a strong rain-shadow effect of the Bernese Alps, the nearby Mediterranean sea influence, and föhn winds on the leeward side, make the Pennine Alps one of the driest regions in the country [54]. Föhn winds blow from the southwest with warm air advected from the Mediterranean or even the Sahara. Having shed its moisture, the dry and warm föhn is then channeled through the northern valleys [65]. Given these physical processes, precipitation within the area varies from <600 mm per year in the lower valleys (~500 m a.s.l.) to more than 2500 mm in very wet mountains [66]. Under such climate conditions, glacier distribution on both sides of the valley considerably differs, starting from 1600 m a.s.l. in the north and 2000 m a.s.l. in the south.

Glaciers play a key role in both the hydrological regime (~10% of the annual runoff) and economic development (e.g., tourism and recreation, hydropower, and irrigation) of Canton Valais. They cover about 570 km² of the landscape (~10% of the catchment) based on the SGI2010, which is about 80% of the total Swiss glacier ice volume [68]. The impact of climate change in terms of glacier evolution and hydrology have extensively been addressed in previous studies [69,70]. According to Fischer et al. [71], the Upper Rhône catchment has lost about 152 km² of the glacierized area since 1973, with a mean geodetic mass balance change of −0.59 m w.e. yr⁻¹ from 1980 to 2010. Thereby, the Rhône will likely be the region of the Alps with the largest land ice remainders (e.g., Aletsch, Gorner and Corbassière (Figure 1)) by the end of the century (~20–30% of the current volume) [68]. Presently, the catchment includes 37 glacier bodies with areas larger than 3 km² that constitute 69% of the total glacierized area and 4.3% of the total number of glaciers in the Upper Rhône catchment. Glaciers smaller than 1 km² account for 90% of the total number of glaciers, while covering 15% of the area. In this study, 849 glaciers within the Bernese and Valais Alps were selected, but we only considered glaciers larger than 0.01 km² for lake delineation. This threshold was applied with respect to the sensitivity of small glaciers to the model parameters.

3. Data

The Swiss Alps has numerous sources of high-quality data, including high-resolution DEMs, historical maps, and a long record of glacier observations. The VOLTA model requires only two input datasets: glacier outlines and a DEM. With the accelerating rate of glacier recession observed in the European Alps during the last few decades, glacier outlines have been consistently updated. In this study, we used two different sources of glacier boundaries and DEMs (see Table 1 and Figure 2). The latest available versions (SGI2010 and a 10-m swissALTI^3D) were used to estimate the glacier ice thickness for current years in order to extract potential lakes, while the second versions (SGI2000 and HDM25 Level 2) were adopted as supplementary data for intercomparison of VOLTA and previous models.

The latest available glacier inventory for the entire Swiss Alps (SGI2010) was derived by Fischer et al. [74]. The accuracy of this inventory was assessed by comparing the extents of clean, snow and/or debris-covered glaciers derived from multiple digitization. The second glacier outlines were extracted from SGI2000, which was created by Paul et al. [73], using Landsat Thematic Mapper (TM) imagery from the years 1998 and 1999. The two versions of DEMs were also obtained from the Swiss Federal Office of Topography (Swisstopo) and used in this study [67,72]. The first DEM, Digital Height Model (DHM25) Level 2, was derived from the Swiss National Map 1:25,000 around 1996. The model corresponds best to the glacier outlines of 2000 and consists of digitized lakes, break lines, and spot heights for interpolation in a symmetric 25-m grid. According to Swisstopo documents [72], the accuracy of this dataset has been determined by comparison of surveyed ground control points (vertical accuracy of 1.5 m for flat areas and 3 m for the Alps). The dataset compounds 1.5 m to 2.2 m mean error of the basis model from the spot heights within the glacierized areas in the Upper Rhône catchment.

The second and more recent DEM is the swissALTI^3D, which represents the surface of Switzerland without vegetation and infrastructure at a 10-m spatial resolution [67]. SwissALTI^3D has been available since 2012 and updated in a 6-year cycle based on aerial photographs. Compared to the previous version (e.g., DHMs), swissALTI^3D was derived by stereo correlation and manual 3D stereo measurements on a 25-cm resolution of 2008–2011 SWISSIMAGE Level 2 aerial ortho-photographs for areas above 2000 m a.s.l. The vertical accuracy of glacier areas in the Upper Rhône catchment at 1σ level is approximately ±1 m to 3 m for the stereo correlation and ±0.1 m to 1 m for manual stereo measurements. Areas below 2000 m a.s.l. were created using airborne laser scanning data with an accuracy of ±0.5 m at 1σ level [67]. In the 2019 edition of swissALTI^3D, the Bernese Alps were tracked with the aerial photos of 2017 and the Pennine Alps were mainly updated by 2016 and partially with 2015 datasets.

Ground Penetrating Radar (GPR) profiles for 42 glaciers were used to assess the accuracy of both modeled ice thickness and basal topography. This dataset was originally compiled as a contribution to the Glacier Thickness Database (GlaThiDa) launched by World Glacier Monitoring Service (WGMS), and is currently available on Glacier Monitoring Switzerland (GLAMOS) [75]. The data were collected in multiple field campaigns from 1999 to 2013, but in this dataset, most of the glaciers were surveyed in the years 2011–2013.

4. Methods

4.1. Glacier Ice Thickness Estimation and Distribution

Linsbauer et al. [50] developed a model based on the shear-stress equation [48] that relates ice thickness to the surface slope. The assumption is that shear stress is constant along the central line of a glacier and the flow is parallel to the bed. Therefore, the ice thickness along this line can be calculated using Equation (1):

(1) $h=\frac{{\tau }_b}{f\rho gsin\alpha }$

where $h$ is ice thickness (m); ${\tau }_b$ is basal shear stress in Pa; $\rho$ is the ice density equal to 900 kg/m³; $g$ is gravitational acceleration (9.81 m/s²) and α refers to the surface slope along the centerline (degree). The shape factor ( $f$ ) refers to the ratio between the width of a glacier and its perimeter, which defines the friction of the glacier with the edge walls. Based on the empirical evidence of alpine glaciers, the values of 0.7 for glacier tongues and 0.9 for wide accumulation zones have been suggested by Maisch and Haeberli [76]. The value of 0.8 is typical for valley glaciers while it can be smaller for other glacier types [48,50,77]. In Equation (1), basal shear stress was derived from:

(2) ${\tau }_b=0.5+159.8\Delta \mathrm{H}-0.435\Delta {\mathrm{H}}^2$

where $\Delta \mathrm{H}$ is the difference in altitude between the highest and lowest point of glacier in kilometers [43]. However, ${\tau }_b$ varies between glaciers, thus a universal value should not be applied, as previous studies show a high uncertainty (up to ±45%) via this method [50]. Driedger and Kennard [78] indicate that glaciers having longitudinal profile >~2600 m reach critical basal shear stress, where the flow cannot adjust the longitudinal profile to a dynamic equilibrium. Therefore, the ice thickness and slope are dependent on that stress. Given this, VOLTA calculates basal shear stress based on a robust elevation band approach:

(3) ${\tau }_b={2.7·10}^4{{\displaystyle \sum }}_{i=1}^n{\left(\frac{A_i}{cos{\alpha }_i}\right)}^{0.106}$

where $A_i$ is the area (m²) of elevation band with 200 m interval, and ${\tau }_b$ is in Pa.

In this study, h was calculated at 100-m successive intervals along the centerlines, mainly to improve computational efficiency and cover the study area properly. The model employs Li et al. [49]’s approach and uses a more physically realistic method to calibrate $f$ according to the local width of glaciers. Following Li et al. [49] and assuming a parabolic cross-section, f can be estimated using the following equation:

(4) $f=1-\frac{1}{1+0.9\frac{w}{h}}$

where w is the half-width of the glacier cross-section at the point of thickness calculation. Equation (1) and Equation (4) can be combined to solve ice thickness including the influence of side-friction (Equation (5)), as shown by Li et al. [49]. It should also be noted that VOLTA (Equation (5)) calculates vertical ice thickness, h, perpendicular to a horizontal x-axis in order to make it more compatible for GIS analysis (i.e., using a tan function as opposed to a sin function):

(5) $h=\frac{0.9w\left(\frac{{\tau }_b}{\rho gtan\alpha }\right)}{0.9w-\left(\frac{{\tau }_b}{\rho gtan\alpha }\right)}$

For further discussion on the methodology for ice thickness estimation (and calculation of the f value) the reader is referred to Li et al. [49] and James and Carrivick [57]. During final processing, the model automatically checks for erroneous perpendicular lengths by the following conditions: (i) if the perpendicular width line intersected another centerline or (ii) if the resulting shape factor, f value (Equation (4)), is < 0.474 (a half-width/centerline thickness ratio < 1 [48]). At points where either of these conditions is correct, ice thickness, h, will be estimated by Equation (1), using the average shape factor (f) along the corresponding centerline. Since ice thickness has sensitivity to surface slope, it may approach to infinity as surface slope tends to zero. This indicates an overestimation of ice thickness for flat regions within the glacier boundary (e.g., tributary confluences) [49,52]. To solve this problem in the model, a minimum slope threshold (α₀) is used, to which any lower slope values are set. In this study α₀ was set at 4°, the default for the VOLTA model [57] and the value proposed by Li et al. [49].

Ultimately, the model interpolates the basal ice thickness to the entire glacier using the Australian National University’s Digital Elevation Model (ANUDEM) 5.3 [79]. Here, the basal points are assumed to be elevation points and glacier outlines are the contour lines where ice thickness values are zero. ANUDEM has the advantage of simulating the concave shape of glacier beds [44,50] and has been frequently referenced as an acceptable interpolation method [49,50,52,57].

4.2. Detection of Glacier Bed Overdeepenings

Figure 2 illustrates the workflow of glacier bed overdeepening delineation. The initial input of the model is an ice-free DEM generated by subtracting distributed ice thickness from ice surface. In theory, overdeepenings can be defined as topographic depressions without a water outlet (Figure 3) and can be abraded and widened by the circular motion of ice, debris materials and subglacial meltwater flow. As a result of this rotational movement, the maximum depth would be the location of maximum erosion from the greatest ice weight.

Given this definition, overdeepenings can be detected easily using GIS hydrology toolset by first identifying the most low-lying cell within a depression with no downslope flow direction [80]. Herein, considering the DEM artifacts and model errors, sinks with depths equal or smaller than the total vertical accuracy of subsurface topography were treated separately. Following that, all depressions were filled to their horizontal pour point levels and overdeepenings were extracted from the difference between the error-filled sinks and the correct-filled sinks on a cell-by-cell basis (Figure 2). If the difference values were less than 0, a true value was returned and assigned accordingly. In the final stage, the depth value of cells within each overdeepening was summarized and the total volume was obtained from the squared cell size multiplied by the sum of the cell depths (defined as the sum of the filled depths minus the sum of margin of errors). Subsequently, the overdeepenings were converted into individual polygons and the summarized values were joined to the corresponding attribute table.

4.3. Classification of Potential Future Lakes

To identify potential future lakes, four morphological criteria, defined by Frey et al. [61] based on ice mechanical behavior and further developed by Magnin et al. [62] and Colonia et al. [81], were investigated in this study area. According to them, a lake development is highly probable where: (i) the surface slope is quite flat (<5°–10°), (ii) a slope break zone is located downslope, (iii) a bedrock threshold and/or decrease in glacier width is evident downstream, and (iv) a crevasse-free zone is followed by heavy crevasses formation indicating different rates of movement (see the details in Figure 10). For each criterion, a value ranging from 1 to 5 was assigned to define its likelihood (Table 2).

4.4. Evaluation of VOLTA and Intercomparison of Ice Thickness Models

We evaluated the VOLTA results in three steps. First, the modeled ice thickness and bedrock depths are compared to the independent GPR measurements provided by GLAMOS for 46 glaciers within the Upper Rhône catchment. The general uncertainty of VOLTA was previously assessed for one of the Swiss glaciers, Unteraar, outside the Upper Rhône catchment [57], but it has not been evaluated for individual glaciers on a point basis.

In the second step, we compared the VOLTA performance with two popular approaches to assess the relative strengths and weaknesses of our model. The ice thickness of the Swiss glaciers was previously measured by numerical flow models using glacier mass balance and ice flow data [52,82,83]. ITEM was first introduced by Farinotti et al. [52] and further refined by Farinotti and colleagues [82]. The new model determines ice thickness based on the glacier mass balance rather than velocity observations in the first model. Despite the difficulty and complexity of input data collection, the evaluation of model with GPR data revealed a high level of accuracy in ice thickness estimation. The second popular method is GlabTop, which is a shear-stress-based model that estimates the ice thickness within 50-m elevation bins based on the Equations (1) and (2). The first version requires manual digitization of branch lines [50], while the second one (GlabTop 2) automatically estimates the ice thickness based on a random selection of DEM cells [63]. The ice thickness values along the centerlines are then interpolated to obtain the distributed ice thickness. This model has been widely used for future lake formation in the European Alps [50,62,84], Himalayas-Karakoram-Tien Shan [63,85,86,87], and Peruvian Andes [81,88]. Here, the depth values of the Swiss glaciers using the above-mentioned approaches were derived from SGI2000 and DHM25 L2. To make these models comparable, we first needed to run the VOLTA using the corresponding datasets. In the final step, the locations, cross-sectional profiles and geometry of VOLTA’s overdeepenings were compared with the GlabTop results.

5. Results

5.1. Ice Thickness and Glacier Volume

In this section, we focus on the ice thicknesses and volumes derived from VOLTA. Figure 4 shows the present ice thickness distribution for all glaciers located in the Upper Rhône catchment, in which ice thicknesses <100 m are dominant.

As reported in Table 3, over 62% of the area (354 km²) is covered by glaciers with ice thickness <50 m, while only ~16% (~94 km²) has ice thickness in the range of 50–100 m. The VOLTA model estimated a total ice volume of 37.17 km³ with a deviation of 12.26 km³ (1σ), which is approximately ±32% of the total ice volume. The deviation from mean is higher for <50 m class, which is largely due to the prevalent small glaciers located in settings with higher elevation and steeper slopes (Table 3). This indicates the sensitivity of these glaciers to the model parameters. The present average ice thickness for the Upper Rhône catchment is ~65 m, and the maximum ice thickness is approximately 578 m. In terms of the glacierized area classes, about 60% of the Rhône ice volume is stored in the 10 largest glaciers (>10 km²), among which Aletsch alone contains ~27% (A ≈ 77 km²) of the total volume. Considering the 27 glaciers >3 km², they contain 84% of ice volume in total.

5.2. Subglacial Topography

The longitudinal profile of the four largest glaciers (Figure 1) is presented in Figure 5, showing the modeled bed elevation along its centerline, which was automatically extracted by VOLTA. The large and thick glaciers have relatively less rugged bed (e.g., Aletsch) than the smaller glaciers. From this sample, the modeled tongue of glaciers lies between 1700 m a.s.l. (Aletsch and Fiescher) to just below 2500 m a.s.l. (Gorner and Corbassière).

Among them, Fiescher’s thin tongue is distinct from the rest of the glaciers. This can be explained by its steep and thin terminus area and relatively gentle-slope and thicker accumulation zone. Unlike Fiescher, Gorner has a thick terminus and a very thin and steep accumulation zone. The other two glaciers have lower ice thickness in higher elevations within the accumulation area, and larger ice thickness in the middle altitudes and their terminus area in the ablation zone (Figure 5). Such differences in ice thickness, depth and bed slope distribution imply the potential locations for future lake formation.

The hypsometric curves of mean ice thickness, volume, and area of glaciers are also represented in Figure 6. The mean ice thickness has a skewed distribution towards the lower elevations (~2100 to 2500 m a.s.l.), whereas the volume and area values are fairly close to a normal distribution with a peak around 3300 m a.s.l. Figure 6 also reveals the increases in volume and mean ice thickness over 2600 m a.s.l., while the glacial area approaches zero at this elevation. This divergence in distributions is likely due to the convergence of several tributaries, which creates a flat area with higher ice thicknesses in large glaciers (i.e., Aletsch), or can be the result of large depressions within the terminus zones.

5.3. Model Performance and Morphological Analysis of Overdeepenings

The overdeepenings beneath the glacier beds were extracted and their corresponding attributes were added based on the methods described in Section 4. The number of predicted overdeepenings varies to a certain extent by the input DEM. To assess the performance of the model, we selected glaciers with a relatively consistent area in both glacier inventories. Here, we found glaciers >5 km² to be more suitable for comparing DEM influence on the modeled overdeepenings. According to Figure 7a,b, VOLTA predicted 125 overdeepenings using swissALTI^3D, and 113 using DHM25 L2. The large differences between the two datasets are particularly evident when depicting the spread of surface area (Figure 7a,b). The largest and deepest overdeepenings are predicted by DHM25 (~280 ha), while both DEMs could detect small overdeepenings in the same range.

More than half of the overdeepenings have a surface area <10 ha, (Figure 7b) and three of them have a surface area above 100 ha according to the both DEMs. In addition, ~25% of the overdeepenings have a maximum depth of less than 15 m and 13 m using the DHM25 L2 and the swissALTI^3D, respectively (Figure 7b). Overall, the obtained overdeepenings represent about 1.8% to 2.1% of the total glacier volume in the Upper Rhône catchment. In relative terms, the swissALTI^3D seems to better predict future lakes due to smaller outliers for the area and mean depth parameters.

The largest overdeepening with an area of 157 ha and a total volume of 96 M m³ is located in the terminus area of Gorner glacier. Given its debris-covered surface, rapid snout degradation and downwasting were observed over the seven past decades (Figure 8). The other two largest glacier overdeepenings are located in the terminus area of Aletsch (A = 137 ha) and Oberaletsch (A = 103 ha). Using both the DEMs, the deepest overdeepenings were found under the accumulation zone of Triftgletscher (D_mean = 70.50 m, A = 36 ha) in Zermatt valley, while the overdeepening of Oberaletsch (D_mean = 70 m, A = 103 ha) holds the most ice volume among the others (~18% of the total volume).

Figure 9 represents the maximum depth of overdeepenings as a function of area and average depth. In general, the larger overdeepenings are expected to be deeper. Our results also confirm this relationship, as the area parameter (Figure 9a) appears to be relatively related to maximum depth (r² ≈ 0.74). Although the maximum depth variability increases for glaciers larger than 0.1 km², the smaller overdeepenings appear to be mainly shallower,. The linear regression between the predicted mean depth and maximum depth (Figure 9b) is also substantially strong (r² ≈ 0.98), with the maximum depth ~2.61 times larger than the mean depth. Key, therefore, is the extent of overdeepenings, which reflects general erosion depths/or power.

5.4. Analysis of Potential Lakes

In the first run (using SGI2010 and the 10-m swissALTI^3D), 171 overdeepenings were predicted and their likelihood of future lake formation was analyzed. Figure 10 represents the modeled overdeepenings for Mont Miné Glacier according to the four determining morphological criteria, described in Section 4.3 (Table 2), highlighting our realistic results. The spatial distribution of overdeepenings shows that overdeepenings primarily develop in certain contexts and locations: (i) cirques; (ii) confluences and valley-constrictions; (iii) specific lithological setting and strength; and (iv) under glacier termini. Among the predicted overdeepenings, a total of 100 have a likelihood of high to very high (A ≈ 872 km², V ≈ 0.56 km³, D_mean ≈ 36 m), while 11 overdeepenings were appeared less likely to be a lake in the future (Table 4, Figure 11). The latter outcomes are the smallest or largest (e.g., Gorner, Aletsch, Oberaletsch) that do not meet the presence of two criteria: bedrock threshold/glacier width change, and crevasses.

5.5. Evaluation of VOLTA and Intercomparison of Ice Thickness Models

The measured bedrock elevation from seven GPR profiles, surveyed in 2003, over Rhône glacier is compared to the predicted bedrock elevation obtained from VOLTA and GlabTop results (Figure 12). The GPR data often show lower elevations, or higher ice thickness, but the VOLTA results are nearly always within the ±24% uncertainty range. GlabTop, by contrast, has substantially underestimated ice thicknesses, as the resulting profiles are located at higher elevations in Figure 12, except at profiles a and g.

The relationship between the measured and predicted bedrock also confirms the less variability in VOLTA’s results (r² ≈ 0.92) versus GlabTop (r² ≈ 0.82) (Figure 12b,c). In particular, at profile d, where the glacier curves, GlabTop fails to accurately delineate the shape of the valley. Although VOLTA models the parabolic shape of glacier beds in good agreement with the shape of GPR profiles, the skewness on the right side of the glacier is a common feature. Comparing the scatterplots of measured and modeled bedrock at GPR measurements, reveals the highest (−97 m) and the lowest average deviation (−19.5 m) at profiles d and e, respectively. The maximum point differences are up to −197 m, mainly located at glacier apices.

Glacier thickness values derived from the other two approaches [50,82] are presented in Table 5 for 11 measured glaciers. The results show that VOLTA has mainly underestimated the glacier ice thickness compared to ITEM, and this underestimation is specifically higher for larger and less steep glaciers. However, VOLTA volumes fall within the range of 95% confidence intervals of ITEM [82] and are nearly close to GlabTop values [50]. This highlights that the VOLTA values are consistent, not too dispersed, and are corresponding well with two significantly different models, i.e., ITEM and GlabTop. Overall, the mean percentage of ice thickness error (MPE) for all measured profiles in the Upper Rhône catchment is around ±22%, of which 28 out of 42 glaciers are underestimated.

5.6. Visual Evaluation of Predicted Overdeepenings

In this section, we present the observations of several glaciers to evidence the robustness of the model for predicting future lakes using the SGI2010 inventory and new swissALTI^3D DEM. These lakes are formed after the generation of glacier inventories, and can provide insights into the model performance for this study area.

The location of future lakes appears to be controlled by the location of stagnant ice development during the glacier retreat. Here, the slope and width along the valley centerline are the two determining factors, where the concave-up form of the valley reaches a region of shallow bed slope, while the steep slopes in the head advance active retreat [89,90]. Moiry glacier in the Valais Alps is a good example to consider (Figure 13), where the rapid frontal retreat over a seven-year time period (2009–2016) has resulted in lake formation. The model was able to capture this overdeepening that stretches over 250 m ahead of the terminus area (Figure 13). The morphological analysis estimates a medium likelihood of lake formation, mainly due to a moderately steep surface slope (~11°) and the absence of crevasses. This overdeepening covers ~4 ha with an average depth of 9 m, continuing to be fed by meltwater input. Figure 13b also shows the upglacier area, where all four morphological conditions are met, and the moulin and wave-ogive features become visible in the upper-bound of forefront crevassing. This upstream predicted overdeepening (A ≈ 6.5 ha, D_mean ≈ 14 m) is located in α < 5°, where the formed crevasses and moulins play a conduit role in focusing the upper surface melt to the head of the developing overdeepening.

Similarly, several small and shallow lakes with a mean depth of ~10 m and a total area of 18 ha have emerged in deglaciated areas since 2010. Figure 14b highlights one of these small lakes (A ≈ 4 ha, D_mean ≈ 7 m) in Lower Theodul Glacier, which used to be one of the Gorner’s tributaries before the frontal retreat. Somewhat surprisingly, this example highlights the possibility of lake formation over steep subglacial beds. Although the model assigns a very low weight to this class of slope, the presence of bedrock threshold/adverse slope is a common characteristic among this group that promotes lake formation.

Figure 15 shows a larger case in Ferpècle Glacier in the Valais Alps where previously several lake outburst floods occurred in 1943 and 1952 [91]. At that time, the glacier tongue was connected to the neighboring Glacier, Mont Miné, but due to the accelerated retreat rate, both glaciers are currently isolated within their valleys. Since the 1950s, the hydrological regime of the valley area has been considerably impacted by the Grande Dixence hydropower. Today, three pumping stations located in the Hérens valley collect meltwater from the Ferpècle and Mont Miné glaciers [92] (Figure 10 and Figure 11). Figure 15b–d demonstrates the temporal evolution of a lake from the collapse of moulins, thanks to the increase of ice movement at the glacier tongue from 2009 to 2016. The lake is presently in the initial stage and the model estimates a coverage of ~10 ha and an average depth of ~27 m.

From a general perspective, the visual comparison between modeled overdeepenings and existing lakes demonstrated the robust performance of our model in the prediction of small to large lakes. No false lakes were found in deglaciated areas, at least in terms of location. An exception is an erroneous area and depth of a lake (A ≈ 9 ha, D_mean ≈ 31 m) located over the steep bedrock of Baltschieder glacier in the Bernese Alps (Figure 16). Although the morphological analysis reveals a low chance of lake formation, likely due to the very steep surface slope (~21°), the break in the slope and presence of bedrock (or adverse slope) have largely contributed to the formation of a small lake.

Close-ups of the site display the oblique direction of the dip that has led to the evacuation of meltwater through two stream channels (Figure 16b,c). This example is a demonstration of tectonic processes (e.g., englacial thrusting) where the tectonic uplifts hinder lake expansion, meaning that the combination of steep slopes and faults provide conduits for fluid flow. The steep and rectilinear morphology of this small subglacial lake also confirms its tectonic origin. Here, the error in the modeled area of this lake can be attributed to the weak performance of the model for steep and small glaciers.

6. Discussion

6.1. Model Performance in Ice Thickness Estimation

The results show that the VOLTA has underestimated the total glacier volumes compared to the ITEM and GlabTop model. The fact is most significant for glaciers smaller than 0.01 km², where the mean ice thickness value for the year 2000 is four times larger in GlabTop results of the entire Swiss glaciers [50]. This can be explained by the sensitivity of this area class to model parameters, where a larger slope distance (αd) reduces the density of ice thickness estimation points, which in turn results in an underestimation of values at the end of each centerline. Therefore, the interpolations approach to zero, especially in terminus zones, and causing high uncertainty for very small-sized glaciers. In general, the volume of small glaciers with areas of extremely high or low slopes is likely to be most sensitive to VOLTA input parameters, as it is evident in the resulting erroneous lakes. On the other hand, the larger glaciers with gentle slope surfaces are the least sensitive. Here, the most effective strategy suggests focusing on more than one approach for different glaciers and regions.

The comparison of ice volume values, returned by VOLTA, with GPR field measurements for different glaciers around the world, shows a range of 25% underestimation and 16.6% overestimation [57]. The results for the Swiss glaciers in the presented case study also confirm a ±22% mean ice thickness error. Considering the ice thickness distribution, VOLTA is in good agreement with field measured profiles as the variations are largely matched, especially where the glacier cross-sections are closer to V-shaped profiles [57]. VOLTA is also capable of capturing relatively accurate maximum ice thickness, mainly close to GPR profiles, for U-shaped valleys. However, the applied ANUDEM algorithm produces V-shaped profiles for corresponding valleys and skews bed topography as a result of centerline deflection, especially at apices. This eventually leads to an underestimation of ice thickness towards the valley walls, yielding less total ice volume for this type of glacier. This model limitation can partly explain the underestimation of VOLTA compared to ITEM (Table 5). For instance, in Corbassière glacier, the derived mean ice thickness value by VOLTA is ~12 m lower than ITEM. However, both VOLTA and GlabTop [50] corroborated that the mean ice thickness values of Corbassière glacier are generally closer to GPR data (MPE ≈ 3%) than ITEM. Overall, the accuracy of VOLTA is highly dependent on the shape of centerlines and valleys, where the average ice thickness values and resulted bed topography are more accurate for straight glaciers with bedrock cross-sections [57].

The Ice Thickness Models Intercomparison eXperiment (ITMIX) project launched by the International Association of Cryospheric Sciences [93], ranks the ITEM and GlabTop fifth and seventh respectively, among the most effective models based on their average performance. According to our evaluation, we expect that VOLTA would be placed somewhere between these two models.

6.2. Overdeepenings and Future Lakes

Typical overdeepenings are prevalent under former termini, where they are enclosed in a moraine context (e.g., Moiry and Gorner). Terminal (or formerly terminal) bedrock overdeepenings are well-documented in the Swiss Alpine Foreland [94]. According to our results, the largest future overdeepening is located in the ablation zone (under the terminus) of Gorner glacier, where both morphologic settings and crevasse patterns allow the development of seasonal lakes. The location of this overdeepening along with more than 40 englacial caves and some marginal tunnels has been confirmed by previous studies [92,95]. The total length of major supraglacial channels (bédières) in this glacier reaches to 22 km. In addition, the observations on Gorner have shown the substantial increase of moulins in terms of number and life cycle since 1985, which could be attributed to either a slower movement or to the ongoing deglaciation phase. Besides the 1980s, our visual evaluation also signifies another turning point for the Swiss glaciers, as most of the lakes started to appear, mainly after the collapse of ice tunnels in the 2010s.

The overdeepenings may also have different origins other than glaciers, as their location can simply explain their origin [24,96,97,98]. Overdeepenings beyond the ELA are often the result of maximum ice velocity, which generally increases at cirques and confluences through ice accumulation and erosion concentration (e.g., Mont Miné). Where lithological settings become weaker, the erosion forces enhance the other conditions for subglacial basin formation. Ultimately, the ice flux decreases in the ablation zone, but the meltwater forces sediment evacuation, and therefore, restores the ice-erosion power. Nevertheless, some are bounded to tectonic structures that progressively deepen by ice-water forcing interactions during successive glacial advances. In general, overdeepenings are considered an equilibrium landform that all glacier beds are eventually favored to reach towards [24]. They develop during glacier retreat, when surface meltwater accesses the bed through crevasses leading to the highest rate in basal sliding, sediment depletion, and erosion power. On the other hand, during the glacier advances when the meltwater decreases, the ice and abundant sediment supplies will rework overdeepenings in pre-existing subglacial basins (e.g., Ref. [99]).

Studies on the morphology of glacial overdeepenings have been limited. Among the few parameters, the relationship between the maximum depth is an important indicator of the erosive power of glacier in the past, and a determining factor in GLOF risk, as the capacity of water storage, in the future. Despite the GlabTop results, from the analysis of the Swiss Alps, the Himalaya-Karakoram region and Peruvian Andes [45], and the empirical approaches suggested by Cook and Quincey [100], VOLTA shows a stronger relationship between the maximum depth and area implying the area-dependent thickness estimations for this study area. Recent observations on the Swiss Alpine foreland and valleys report a strong relationship between the lithology of glacier bed and the geometry of overdeepenings, as the larger, wider, and shallower overdeepenings tend to be dominant in very low erosional-resistant bedrock [101]. However, the modern overdeepenings in the Swiss Alps, especially in the Valais region, are predominately located in medium to highly resistant bedrocks. The findings further highlight the headward erosion and the efficiency of sediment evacuation, which all depends on the amount of available meltwater during the melt season.

6.3. Limitation of Methodology and Future Applications

The resulting bed topography highly relies on the accuracy of the input DEMs. In general, high-quality data and DEM resolution allow an accurate calculation of surface slope, which ultimately leads to the robust locating of overdeepenings. Frey and Paul [37] confirm that a coarser DEM resolution decreases maximum and increases minimum elevation since these two parameters depend on the individual cell values which are strongly influenced by the resolution of DEM. Herein, the comparison of 25-m and 10-m DEMs shows a considerable influence of DEM resolution on the area of overdeepenings. We also compared the VOLTA’s overdeepenings with those extracted by [62], using GlabTop 2 [63], and 20 m to 25 m resolution of IGN and ASTER DEMs, respectively. In contrast with VOLTA, not only the location and area but also the depth of overdeepenings varied significantly by the choice of DEM. This can be induced substantially by the artifacts present in the ASTER DEM resulting in less variability of ice thickness and thus shallower overdeepenings, especially in flat regions. Thus, this supports the impact of input parameters on the general model performance. Further, Frey and Paul [37] recommend that the acquisition date of DEM should not be older than the acquisition date of glacier outlines, especially for the glaciers that are retreating fast. We refer the partial uncertainty present in our results to the uncertainties of the glacier outlines, the time of the DEMs and field data acquisitions, especially for the possible inconsistency of SGI2000 and DHM25, which were used for the evaluation of the model performance.

The application of four morphometric criteria on historical topographic maps for lake prediction has previously shown a 75% possibility of lake formation, where all criteria are present [7,45]. However, there are other controlling factors in lake formation. Ongoing global warming and change in precipitation patterns have caused the rise of snow line altitudes exposing the clear-ice surface. In the context of prolonged ablation season, the glacier melt has accelerated as well, resulting in feedback mechanism change, e.g., albedo [102,103]. Recently, the darkening glacier surface has been observed and discussed, especially for the Swiss glaciers. According to Naegeli and co-authors [104,105], two glaciers, Brenay and Ferpècle (Figure 11 and Figure 15b) in the Valais Alps, have shown significantly positive trends on average glacier-wide bare-ice albedo since 1999. In addition, the presence of thick debris over the glacier surface can reduce the rate of ablation and delay the time of lake formation if the glacier is not calving into water [106]. Thus, there is a possibility that the predicted overdeepenings, mainly located in glacier tongues and covered by thick debris, transform into outwash plains or shallow lakes/depressions while being filled with debris from the surrounding margin valleys.

Nevertheless, our knowledge of overdeepening and lake formation modeling is limited. To address such shortcomings, initiatives involving repeated field measurements, data collection on the hydrological configuration of overdeepenings, bathymetric study of current lakes, and reconstruction of former glaciers for integration into the modeling process would be the opportunities for future research. Currently, the GPR points are inconsistently distributed over the glacier bodies. They are mainly located in termini, while lacking over the steep parts, and crevasse zones of glacier that play a key role in lake development. Furthermore, direct bed measurements are subject to inaccuracies during the data collection and subsequent data manipulation process. The major source of errors in estimation of bedrock depth using the helicopter-borne technique is due to an inaccurate estimation of GPR wave velocity and picking errors of the surface and bedrock reflections [107]. The variation in helicopter speed and estimation of antennae altitude from the reflected signals of glacier surface during the data collection can lead to spatial irregularities and false bedrock positioning, especially in flat and thick areas of the glacier (i.e., confluences). Due to such challenges, a direct validation is not fully possible and requires further efforts, such as the Glacier Ice Thickness Database (GlaThiDa) [108,109], for the future generation of ice thickness modeling and subglacial basin studies.

7. Conclusions

The increasing number and volume of lakes in glacier forelands have raised the necessity of glacier bed modeling. This study investigated the favorable conditions and sites for future lake development within the temperate glaciers of the Swiss Alps. In addition, the efficiency of VOLTA for the prediction of subglacial topography was examined by intercomparing the results with field measurement data and previous models for this region.

The model estimated a total ice volume of 37 ± 12 km³ (1σ) and an average ice thickness of ~65 m for the total glacierized area (~577 km²) using the high-quality inventories. Although VOLTA can estimate maximum ice thickness for U-shaped valleys close to the field measured profiles, the modeled V-shaped profiles, using the ANUDEM algorithm, for this form of valleys underestimate the ice thickness towards the valley walls. As a result, the accuracy of VOLTA is highly dependent on the shape of the valley and centerlines. This fact implies that the average error in ice thickness (~±22%) is not completely related to the model robustness but also varies from one glacier to another. Comparing the results with ITEM and GlabTop models highlights the tendency of the model in underestimation of average ice thickness, especially for small glaciers located in higher elevations over steep bedrock. However, VOLTA uses an enhanced centerline generation algorithm by including a side drag parameter in the equation, which reduces the uncertainty of glacier bed elevation and therefore results in more accurate cross-sectional profiles than GlabTop.

In this study, the model could identify a total of 171 overdeepenings, in which 100 of them were classified as highly likely to be a glacial lake in the future. The total volume of high potential class is 2.39 km³ with an average depth of 22 m. The largest overdeepening in the Upper Rhône catchment is located in the terminus zone of Gorner glacier with a surface area of 157 ha. We also observed and reported several examples of lakes formed in deglaciated areas between 2010 and 2019, which confirms the robustness of our model in lake/overdeepenings prediction. Such observations provide a compelling case for modifying and applying the mentioned models to understand the glacial evolution at regional scales.

Overall, bed overdeepenings are able to change the entire dynamic behavior of glaciers. Foremost, with ongoing anthropogenic global warming, the presence of overdeepenings and their corresponding processes should be taken into account for modeling the glacier landscapes. However, the studies on the implications of such processes are still lacking, which makes the modeling process difficult. The approaches like VOLTA may predict the location of lakes accurately to some extent, but volume and area still face levels of uncertainty, largely for small steep glaciers. Future research requires addressing these uncertainties in the model by focusing on developing algorithms capable of simulating valley shapes more precisely, in order to avoid the underestimation of ice thickness adjacent to steep valley sides.

To summarize, the model presented here can be an alternative to current models which use sparsely available field-based parameters such as mass balance measurements [85,110], or surface velocity [111], and models that require manual extraction of centerlines [50].

Author Contributions

Conceptualization, S.G., J.L.R.J.; methodology, S.G., W.H.M.J.; software, S.G.; validation, S.G., A.B., L.S.; formal analysis, S.G., W.H.M.J.; investigation, S.G.; resources, S.G., A.B., W.H.M.J.; data curation, S.G.; writing—original draft preparation, S.G.; writing—review and editing, J.L.R.J., T.J.B., D.R.B., W.H.M.J., A.B., S.G.,and L.S.; visualization, S.G.; supervision, J.L.R.J., T.J.B., and D.R.B.; funding acquisition, S.G., J.L.R.J., A.B., and L.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research is funded by Texas State University through a Doctoral Research Support Fellowship.

Acknowledgments

We thank Paul, Fischer and Linsbauer for sharing SGI2000 and Rhône Glacier data. Magnin and Bauder for sharing their research results, Weaver and Klimeš for early feedback, and Chenoweth for final revision. The Swiss glacier inventories and ice thickness measurements are available on GLAMOS and the DEMs and images can be obtained from Swisstopo.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Overview map of the Upper Rhône catchment in Canton Valais of Switzerland. Purple outlines show the extent of glaciers derived from SGI2010, lying over a hillshade of swissALTI3D Digital Elevation Model (DEM) [67]. The colored glaciers were used for the validation of bed topography along the extracted centerlines by Volume and Topography Automation (VOLTA) in Figure 5: Great Aletsch (ALE); Fiescher (FIE); Gorner (GOR); Corbassière (COR). The inset map features the configuration of the Upper Rhône catchment at a national scale and the abbreviations outside the Swiss (CH) borders refer to the neighboring countries: Italy (IT); and France (FR).

Figure 2. Schematic illustration of the workflow and processing steps of the proposed VOLTA approach for ice thickness estimation (a), and the geoprocessing approach for extracting overdeepenings (b).

Figure 3. Two-dimensional cross-section of an overdeepening (a) and the process of lake formation after glacier retreat (b). The blue arrow in (a) shows the direction of ice flow.

Figure 4. Ice thickness distribution of the Upper Rhône catchment glaciers (Data source: SGI2010, 10-m swissALTI3D). A statistic summary of ice thickness classes is presented in Table 3 The labeled glaciers are additionally referenced in Table 5.

Figure 5. Longitudinal profile of modeled beds (dashed line) of four largest glaciers in the Upper Rhône catchment (Data source: SGI2010, 10-m swissALTI3D). The solid line shows the present glacier surface.

Figure 6. The hypsometric distribution of glacier mean ice thickness, area and volume, averaged over each 100-m elevation interval (Data source: SGI2010, 10-m swissALTI3D).

Figure 7. Effect of DEMs (10-m swissALTI3D and 25-m DHM25 L2) on morphometric parameters and overall model performance for overdeepenings beneath glaciers >5 km2; and summary statistics of the corresponding overdeepenings presented in logarithmic (10 (a) and 2 (b)) scale for better display.

Figure 8. Gorner glacier tongue and the largest predicted overdeepening (dark blue polygon). Red arrows highlight the snout change since 1946 and the dashed ellipses emphasize the temporal terminus degradation and formation of an ice cave/tunnel (Data source: Satellite Image © CNES/SPOT Image/Swisstopo, NPOC (a,b), and Google Earth (c,d)).

Figure 9. Scatter plots of maximum depth as a function of area (a) and average depth (b) presented in semi-log (10) scale for better visualization. Intercepts are set to 0.

Figure 10. Morphological criteria that determine the likelihood of future lakes.The modeled overdeepenings (blue polygons) beneath the accumulation zone of Mont Miné Glacier located in the Hérens valley of the Valais Alps are presented (a) (Data source: VOLTA and SGI2010). These cirque-type overdeepenings are separated by distinct valley steps. The red arrow shows the location of a GLOF and the red circle marks the Ferpècle pumping station and water intake; (b) a schematic diagram (yellow dashed outline) illustrating a potential lake site close to equilibrium line altitude (ELA) where slope increase along the glacier flow direction (black arrow showing the angle change), decrease in glacier width (red arrow), and transition from crevasse-free area (red circle) to heavy-crevasse formation (blue squiggles and red ellipse) are clearly visible.

Figure 11. The spatial distribution of modeled overdeepenings in the context of hydropower development in the Bernese Alps (a) and the Valais Alps (b). The dashed blue polygons show the location of the largest overdeepenings. Great Aletsch (ALE); Ferpècle (FRP); Gorner (GOR); Moiry (MRY); Brenay (BRN); Oberaletsch (OAL); Baltschieder (BAL); Mont Miné (MMN).

Figure 12. Comparison of cross-sectional bed profiles (Pa to Pg) at Rhône glacier derived from VOLTA (blue dots), GlabTop (yellow dots) and field measurements (orange dots). The map shows the locations of the GPR measurement points and modeled overdeepenings using VOLTA and GlabTop (a). Scatter plots (b) and (c) evaluate the relationship between the modeled and measured bedrock values. (Coordinate System: CH1903 LV03; Survey Date: 04/26/2003 [82]).

Figure 13. Temporal evolution of a lake following Moiry Glacier recession in the Valais Alps (a–c). The red solid arrow marks the location of glacier snout in 1946. (a), which is currently ice-free (b–e). The longitudinal and frontal views of the Moiry Valley are demonstrated in (d–e). The notched arrows highlight the location of two moulins over the glacier surface. The blue polygons are the predicted overdeepenings and the gray polygon in (e) shows the glacier outline in 2010. (Data source: Satellite Image © CNES/SPOT Image/Swisstopo, NPOC (a), and Google Earth (b–e)).

Figure 13. Temporal evolution of a lake following Moiry Glacier recession in the Valais Alps (a–c). The red solid arrow marks the location of glacier snout in 1946. (a), which is currently ice-free (b–e). The longitudinal and frontal views of the Moiry Valley are demonstrated in (d–e). The notched arrows highlight the location of two moulins over the glacier surface. The blue polygons are the predicted overdeepenings and the gray polygon in (e) shows the glacier outline in 2010. (Data source: Satellite Image © CNES/SPOT Image/Swisstopo, NPOC (a), and Google Earth (b–e)).

Figure 14. Overdeepenings (blue polygons) at the tongues of Gorner, Lower/Upper Theodul Glaciers in the Valais Alps (a). Red circle shows the presence of a small-sized lake that was not detected by the model and red dashed square and yellow arrow mark the location of emerging lakes (b,c) within the predicted overdeepenings at Theodul Glaciers.

Figure 15. Temporal evolution of a potential lake (blue polygon) following the collapse of moulins in 2016 and the rapid recession of tongue since 1983 (red solid arrow (a)) at Ferpècle Glacier in the Valais Alps. The lake is at initial stage by two sinking moulins (red notched arrows in (c)). (Data source: Satellite Image © CNES/SPOT Image/Swisstopo, NPOC (a,c), and Google Earth (b,d)).

Figure 16. False estimation at Baltschieder glacier in the Bernese Alps (a), including a close-up of a small bedrock-dammed lake (b), and the two outlets draining the meltwater (c) (Data source: Satellite Image © CNES/SPOT Image/Swisstopo, NPOC).

Table 1

Details of acquired DEMs and glacier inventories.

Attributes	DEMs
	DHM25 Level 2	swissALTI3D
Resolution	25 m	10 m
Date of acquisition	1996–2000/2001	2015–2016–2017
Vertical accuracy	1 m–3 m	0.1 m–3 m
Data source	Swisstopo [72]	Swisstopo [67]
	Glacier outlines (SGIs)
	SGI2000	SGI2010
Data source	Landsat TM4-5	Orthophotos
Date of acquisition	1998–1999	2008–2011
Area	>5 km2	All
Corresponding DEM	DHM25 L2	swissALTI3D
Reference	Paul et al. [73]	Fischer et al. [74]

Table 2

Four morphological criteria for determining the likelihood of a lake.

Criterion	Low	Very Low	Medium	High	Very High
Slope	≥20°(>20%)	15°–20°(>50%)	10°–15°(>50%)	5°–10°(>50%)	<5°(>50%)
Break in slope	≤5°(>20%)	5°–10°(> 50%)	10°–15°(>50%)	15°–20°(>50%)	>20°(>50%)
Bedrock/Glacier width change	Absent				Present
No crevasse → heavy crevasse	Absent				Present

Table 3

Ice thickness and Area classes, including Mean Slope (α), Maximum Ice Thickness (h_max), and Mean Ice Thickness (h_mean) (Data source: SGI2010, 10-m swissALTI^3D).

Ice Thickness (m)	Area (km2)	%	Volume (km3)	%
0–50	354.20	62.20	7.70 ± 5.31	20.70
50–100	94.50	16.60	6.80 ± 1.40	18.20
100–200	81.30	14.30	11.32 ± 2.21	30.40
>200	37.80	6.50	11.35 ± 3.33	30.50
Total	577	100	37.17 ± 12.26	100
Area (km2)	∑ Area (km2)	α (°)mean	hmax (m)	hmean (m)
<0.01	0.98	39	15	13
0.01–1	85	37	43	34
1–3	85	26	311	35
3–10	156	21	333	56
>10	241	20	578	94

Table 4

Modeled future lakes and their corresponding likelihood based on the four morphological criteria. The total area, volume and average depth are also calculated for each class.

Likelihood	No. Potential Lakes	Area (km2)	Volume (km3)	Depthmean (m)
Very low	11	5.37	0.028	17.10
Low	16	69.94	0.029	12.96
Medium	47	325.84	0.187	17.00
High	49	470.87	0.353	18.97
Very high	51	401.87	0.207	17.54

Table 5

Glacier parameters including Area (km²), Shear-stress (τ_b), Mean Slope (α), for glaciers > 5 km² with direct ice thickness measurements, and comparison of the individual glacier Maximum Ice Thickness (h_max) and Mean Ice Thickness (h_mean) derived from the VOLTA model (DHM25) with ITEM and GlabTop (using two methods of interpolation: Inverse Distance Weighting (IDW) and Topo to Raster (TTR)) for the SGI2000 outlines.

Glacier *	Area (km2)	τb (kPa)	α (°)	Aspect	hmax (m)	hmean (m)
						VOLTA	ITEM	GlabTop (IDW/TTR)
Aletsch (ALE)	82.15	185.04	14.60	SE	584.13	133	187	140/160
Gorner (GOR)	55.80	172.19	20.46	NW	468.30	72	106	87/96
Corbassière (COR)	16.70	155.54	19.40	NW	290.70	80	92	97/110
Rhône (RHO)	15.76	153.81	15.19	SW	370.55	104	132	99/112
Zinal (ZIN)	14.02	153.49	24.00	NE	219.90	50	66	56/64
Allalin (ALL)	9.27	147.42	18.89	NE	260.58	50	91	62/74
Mont Miné (MMN)	8.98	147.31	15.90	W	241.27	74	76	79/88
Lang (LNG)	9.36	145.51	21.64	S	214.47	52	60	67/71
Mont Durand (MDR)	5.94	140.65	20.34	E	224.67	72	59	67/73
Trient (TRN)	5.99	140.11	17.74	NW	232.20	80	67	76/86
Giétro (GIE)	5.42	138.41	14.50	NW	270.29	104	88	99/104

* All glaciers have direct ice thickness measurements.