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

: 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 3 (1 σ ) was estimated for the Upper Rh ô ne catchment. The comparison of VOLTA performance indicates a stronger relationship between measured and predicted bedrock, conﬁrming the less variability in VOLTA’s results (r 2 ≈ 0.92) than Glacier Bed Topography (GlabTop) (r 2 ≈ 0.82). Overall, the mean percentage of ice thickness error for all measured proﬁles 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.


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 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.

Study Area
The Upper Rhône catchment is located in Canton Valais in southwestern Switzerland (Figure 1), covering a total area of ~5290 km 2 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.  [67]. The colored glaciers were used for the validation of bed topography along the extracted centerlines by Volume and Topography Automation (VOLTA) in  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 2 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 2 of the glacierized area since 1973, with a mean geodetic mass balance change of −0.59 m w.e. yr −1 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 2 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 2 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 2 for lake delineation. This threshold was applied with respect to the sensitivity of small glaciers to the model parameters.

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. Paul et al. [73] Fischer et al. [74] 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, Remote Sens. 2020, 12, 3443 5 of 28 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 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 Level2 aerial orthophotographs 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.

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): 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.

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): where h is ice thickness (m); τ b is basal shear stress in Pa; ρ is the ice density equal to 900 kg/m 3 ; g is gravitational acceleration (9.81 m/s 2 ) 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: where ∆H is the difference in altitude between the highest and lowest point of glacier in kilometers [43]. However, τ 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: where A i is the area (m 2 ) of elevation band with 200 m interval, and τ 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: 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): 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 (α 0 ) is used, to which any lower slope values are set. In this study α 0 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]. 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.

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).
No crevasse→ heavy crevasse Absent Present

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.

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.
Remote Sens. 2020, 11, x FOR PEER REVIEW 9 of 28  Table 3 The labeled glaciers are additionally referenced in Table 5.

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  Table 3 The labeled glaciers are additionally referenced in Table 5.
As reported in Table 3, over 62% of the area (354 km 2 ) is covered by glaciers with ice thickness <50 m, while only~16% (~94 km 2 ) has ice thickness in the range of 50-100 m. The VOLTA model estimated a total ice volume of 37.17 km 3 with a deviation of 12.26 km 3 (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 2 ), among which Aletsch alone contains~27% (A ≈ 77 km 2 ) of the total volume. Considering the 27 glaciers >3 km 2 , they contain 84% of ice volume in total.

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. 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 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.
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.

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 2 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 3 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 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.   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 2 ≈ 0.74). Although the maximum depth variability increases for glaciers larger than 0.1 km 2 , 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 2 ≈ 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.  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 2 ≈ 0.74). Although the maximum depth variability increases for glaciers larger than 0.1 km 2 , 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 2 ≈ 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.

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 2 , V ≈ 0.56 km 3 , 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.  Table 4). 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.

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 2 ≈ 0.92) versus GlabTop (r 2 ≈ 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.

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.
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.

Visual Evaluation of Predicted Overdeepenings
In this section, we present the observations of several glaciers in Figures 13-16, evidencing 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.
Remote Sens. 2020, 11, x FOR PEER REVIEW 17 of 28 Figure 13. Temporal evolution of a lake following Moiry Glacier recession in the Valais Alps (a-c).   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 13. Temporal evolution of a lake following Moiry Glacier recession in the Valais Alps (a-c).      [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] (Figures 10 and 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.  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   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

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 2 , 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.

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.

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 (Figures 11 and 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.

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 3 (1σ) and an average ice thickness of~65 m for the total glacierized area (~577 km 2 ) 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 3 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].