Integrating GPR and ice-thickness models for improved bedrock detection: the case study of Rutor temperate glacier

Estimating the ice volume contained in glaciers is a topic of increasing interest, because the cryosphere has rapidly evolved during the last decades of global warming. Many disciplines collaborate to study the global warming impacts on glaciers. From the perspective of a geophysical research team, we examine the advantages of integrating glaciological ice-thickness models into the workflow for geophysical data processing.

In temperate glaciers, the widespread englacial water content often challenges the analysis of Ground Penetrating Radar (GPR) data, the most commonly used geophysical technique for measuring ice thickness. In fact, recognizing the ice-bedrock interface is hard where the englacial water content generate high levels of scattering in the GPR data. A past GPR survey estimated that the Rutor glacier (European Alps) stored about 150 million m³ of ice in 2008. However, this estimate proved unrealistic after analyzing the geodetic mass balance of the following decade. Therefore, we analyzed new GPR measurements on the same glacier, which highlight the difficulty in identifying the true bed reflection. On these data, we tested the idea that ice-thickness models can help the GPR data analyst to better recognize the ice-bedrock interface in the radargrams. We selected four models, OGGM, GlabTop2-Py, Original-GlabTop2 and GlaTE, that estimate bedrock topography starting from surface topography, following principles of ice flow theory, ice dynamics and mass conservation. Combined visualization of the GPR and model data in 2D and 3D helped the analyst to manually select the ice-bedrock interface. This proved useful especially where the GPR data was scattered and interpretation was uncertain. The GPR data were then used to constrain one of the models, GlaTE, to produce an ice-thickness map that is the result of both model estimates and GPR information. Without this modelling context, the GPR data could be misinterpreted, and the resulting GPR-derived ice-thickness estimates might then be used to constrain a subsequent model in a way that would introduce significant errors. According to this methodology, the Rutor glacier stored about 450 million m³ of ice in 2021, about three times the previous estimate. The workflow is openly available in the section “Code and data availability” and may improve future GPR surveys of temperate glaciers, especially when facing scattered data due to englacial water content or other sparse reflectors such as debris. More accurate ice-thickness estimates will improve local studies and provide better calibration data for regional studies.

1 Introduction

Temperate valley glaciers are a characteristic feature of mountain chains in temperate zones, such as the European Alps, and have a pivotal role in the hydrological cycle (Milner et al., 2017). Studying the glacier processes and changes linked to the current global warming often requires a reliable estimate of their bedrock topography, crevasses, cavities, debris, and englacial water (Haeberli et al., 2019). Temperate glaciers are characterized by containing temperate (or “warm”) ice at the melting point. Indeed, colder ice zones often exist in temperate glaciers, but they are limited to specific areas or seasons, e.g., at the highest altitudes, or where the ice is thin and cools down quiclky during the cold season (Suter et al., 2001; Reinardy et al., 2019). Temperate glaciers at the pressure melting point are primed for rapid meltwater production upon small energy or heat inputs.

A high englacial water content in glaciers can challenge the interpretation of bedrock returns from Ground Penetrating Radar (GPR), one of the most used techniques for surveying ice masses (Colucci et al., 2015; Forte et al., 2015; Urbini et al., 2019). Since englacial water has very different dielectric permittivity compared to ice, it reflects and scatters the electromagnetic wave, hindering the signal from traveling to the bedrock and then being detected on its way back (Reinardy et al., 2019). Smaller-scale heterogeneities, such as small fractures or sediment grains smaller than about a quarter wavelength, generate weak or undetectable responses, but their presence influences the signals as they pass by. The heterogeneities extract energy and scatter the signal in all directions (Jol, 2009). Langhammer et al. (2019b) studied the effects of antenna orientation on detection of the bedrock reflection. Challenges in detecting basal returns over temperate glaciers have been studied by Rutishauser et al. (2016), who analyzed a large set of GPR data acquired on Swiss glaciers and found that depending on the specific glacier, the bedrock interface could only be successfully detected in 12 %–69 % of the data due to scattering. GPR signal scattering rarely occurs in arctic cold glaciers, but when it is detected in some areas, it may be evidence of temperate conditions and englacial water (Karušs et al., 2022). Another common source of noise in GPR sections is the presence of debris at the glacier surface, as noted by Colombero et al. (2019), or at depth, due to phenomena such as adfreezing and entrainment of sediments into the basal ice layer (Weertman, 1961). Air bubbles trapped in ice cause scattering of GPR signal, which helps differentiate between various types of ice, such as firn and superimposed ice (Langley et al., 2009).

To study and address the issue of detecting bed returns obscured by scattering noise in temperate glaciers, we focus on the Rutor Glacier (RGI ID: RGI60-11.03039), the third-largest glacier in the Aosta Valley, Southern Europe Alps. The Rutor massif develops on multiple terraces, over which the Rutor Glacier has advanced or retreated during the last centuries because of the climate variations, shaping the basin morphology, as reviewed by Vergnano et al. (2023). Nowadays, the Rutor Glacier covers an area of 7.5 km² and it is situated in the municipality of La Thuile. It has retreated, since about 1860 (when it was at the Little Ice Age maximum extent), to an upper terrace, forming new proglacial lakes. Its tongues are situated on the northern side of the glacier at about 2550 m a.s.l. (the eastern tongue) and 2650 m (the central and the western tongues). Its highest elevation is about 3440 m, near the southern margin.

In recent decades, the glacier has experienced an accelerated retreat, linked to climate warming (Corte et al., 2024; Gizzi et al., 2022). By differencing 2021 and 2008 DEMs, pixel by pixel, we observed the changes in the glacier surface topography from 2008 to 2021. This simple analysis showed a loss of more than 20 vertical meters of ice in about $\mathrm{1}/\mathrm{4}$ of the current glacier area (especially in the tongue area) for a total ice volume loss of about 100 million m³. The 2008 DEM was retrieved from the regional cartography (Val d'Aosta Region, 2008); the 2021 DEM is available in an open repository described by Corte et al. (2024). Previous GPR surveys (1996–2006) reported an average thickness of the glacier of 17.5 m (a volume of 150 million m³ over an area of 8.5 km²) (Villa et al., 2008, 2007). However, this estimate does not agree with the ice lost from 2008 to 2021. If the previous ice volume estimate was correct, Rutor Glacier would have lost $\mathrm{2}/\mathrm{3}$ of its volume in just one decade, meaning that only 50 million m³ are left. However, while shrinking is evident, the remaining ice volume is likely larger given the relatively small area reduction in that decade (see Fig. 2). We hypothesize that the ice volume of 150 million m³ is underestimated due to problems in interpreting scattered GPR data. The analysis performed on the same glacier by Viani et al. (2020) presents similar considerations about the challenge of interpreting the correct bedrock geometry. Moreover, many examples of other alpine glaciers with similar areas consistently show greater ice thickness (Grab et al., 2021).

In this study, the Rutor Glacier is investigated with two new GPR datasets, acquired in May 2012 with a helicopter-based survey (Morra di Cella, 2024) and in May 2022 with a ground-based survey. These new datasets reveal high scattering of the radar signal over most parts of the glacier, demonstrating the difficulty in detecting the ice-bedrock interface. The scattering zone is often located at a depth of around 20 m, and may easily be misinterpreted as the ice-bedrock interface, possibly explaining the previous doubtful estimates (Villa et al., 2008).

To address this problem encountered on the Rutor Glacier, but common to other temperate glaciers, this study explores the idea that ice-thickness models may aid in the interpretation of scattered GPR sections, especially to fill the gaps in the most scattered survey sections. Those algorithms require, as input, the glacier surface topography. Some models have obtained much acknowledgment in recent years and were reviewed during the ITMIX project (Farinotti et al., 2017), in which many glacier models were compared on the same set of glaciers. One of the conclusions of that research is that the average of different models is a more reliable ice-thickness estimate than a single, finely-tuned model. Such models have been used to perform regional studies, for example, on the Austrian glaciers (Helfricht et al., 2019) and on Swiss glaciers (Grab et al., 2021), the latter including large-scale GPR surveys.

The four ice-thickness models employed in this work are: OGGM (Maussion et al., 2019), GlaTE (Langhammer et al., 2019a), the original GlabTop2 (Frey et al., 2014), and its open-source implementation Glabtop2-Py. The ice thickness is predicted by using the four models. During the manual selection of the ice-bedrock interface in the GPR data, the results from the models are superimposed on the radargram to help identify the most likely ice-bedrock interface, often submerged by noise due to the englacial water content. The GPR data are then used to constrain a second run of the GlaTE model, to provide a final model of the bedrock topography.

2 Methods

We present here an overview of the workflow, enumerating the main steps of the proposed methodology. Then, the next paragraphs analyze the methodology in more detail. For reasons of space, we did not detail here all the technical steps needed to reproduce the workflow. The interested reader is encouraged to download the supplementary materials cited in the section “Code and data availability” (hereafter referred to as Supplement), organized with the same logic used here, which contain “readme.txt” files providing technical help.

1. Run four models (Original-GlabTop2, Glabtop2-Py, GlaTE, OGGM), which, based on the surface topography, estimate the ice thickness. We calculated the ice-thickness average and standard deviation of the four models, to have an overview of the ice-thickness uncertainty.

2. Collect and analyze new GPR datasets from two acquisitions, a helicopter-based survey made in 2012 and a ground-based survey made in 2022. The GPR data is often scattered and difficult to interpret.

3. Extract the ice thickness estimated by the models in the same locations of the GPR paths, thanks to the v.sample tool of QGIS, GRASS plugin (QGIS Development Team, 2021), using a bilinear interpolation method.

4. Visualize the GPR data and the model estimations together, in 2D with RGPR, (Huber and Hans, 2018) and 3D with Paraview (Ayachit, 2015).

5. Manually select (“pick”) the bedrock reflector in the GPR data with the help of the estimations from the three models, to reduce the chance of misinterpreting bedrock reflections. For analysis and discussion purposes, we divided the pickings into “sure” pickings, “guided by model” pickings, and “probably incorrect” pickings.

6. Run the GlaTE model for a second time, constraining the estimations with the GPR data.

7. Run the GlaTE model multiple times, varying the input glaciological parameters, to perform a simple sensitivity analysis on both the unconstrained and the constrained-by-GPR version of the model.

8. Compare the results with literature products.

2.1 Ice thickness models

The GlabTop2, GlaTE and OGGM models used in this study require, as input, a digital elevation model (DEM) of the glacier surface, and they estimate the ice-thickness based on ice flow theory, ice dynamics and mass conservation. In this study, these models were employed for two objectives. First, to provide an initial estimate of the ice thickness distribution of the Rutor Glacier, in order to help identify the ice-bedrock interface in the GPR sections. Second, the GlaTE model was employed to provide a final estimate of the Rutor Glacier ice thickness, based on both model estimates and GPR data constraints.

The models required additional input parameters (e.g. ice density, Glen's exponent, etc.). These were checked for consistency with the Rutor Glacier study area, but unless otherwise stated, the default values from similar alpine glacier studies by the model developers were used. This choice could introduce some errors: these parameters are not transferable between glaciers (Zekollari et al., 2022). For this reason, we performed a sensitivity analysis on the final model (see the specific paragraph).

The DEM and the glacier margin, manually drawn based on a high-resolution orthophoto, are retrieved from a topographical survey carried out in 2021 (Macelloni et al., 2022; Corte et al., 2024). With the warp/reproject tool of QGIS software, using a bilinear-based triangulation method, the DEM was undersampled to 20 m resolution, for computation-time optimization and because using a too-fine resolution was found to produce misleading structures in the resulting bedrock topography.

The features of the three models are described below.

2.1.1 GlabTop2

The GlabTop2 (Glacier bed Topography 2) model assesses the distribution of ice thickness in a glacier starting from a DEM file and a mask file containing the margin of the glacier (Frey et al., 2014). It employs an algorithm first developed by Linsbauer et al. (2012), but slightly modified to avoid the laborious process of manually drawing branch lines:

$$\begin{array}{}\text{(1)}& h_{\mathrm{f}}={\displaystyle \frac{\mathit{\tau }}{f\cdot \mathit{\rho }\cdot g\cdot \sin \left(\mathit{\alpha }\right)}}\end{array}$$

where h_f is the mean ice thickness along the central glacier flow line, f is a shape factor, fixed at 0.8 (according to Haeberli and Hoelzle, 1995; Frey et al., 2014), τ the shear stress at the glacier base, calculated with an empirical formula based on the elevation range of the glacier ΔH (Eq. 2), and expressed in kPa, ρ the ice density (900 kg m⁻³) (Langhammer et al., 2019a), g is the gravity acceleration, and α is the surface topography slope:

$$\begin{array}{}\text{(2)}& \mathit{\tau }=\mathrm{0.5}+\mathrm{159.81}\mathrm{\Delta }H-\mathrm{43.5}\mathrm{\Delta }{H}^{\mathrm{2}}\end{array}$$

The first processing step employed by the algorithm is an initial approximation of ice thickness in some random cells of the domain, based on the surface slope of a sufficiently large buffer zone around the cell. Then, the ice thickness of the remaining cells is estimated with a simple inverse distance weighting algorithm. The two steps are repeated for n times, and the results are averaged.

As the original GlabTop2 code is written in a closed-source environment, we also used an implementation of GlabTop2 in Python, available open-source in a public repository (Terink, 2018). Further details are provided in the appendix of Frey et al. (2014). Glabtop2-Py is a partial implementation of the original code, and we noted a poorer performance near the glacier margins. In the Supplement, the GlabTop2-Py model code adapted for this case study is available to ensure full reproducibility. The config.cfg file contains the input parameters used in this study.

2.1.2 GlaTE

GlaTE (Glacier Thickness Estimation) is also based on Eq. (1), but with a different estimate of the shear stress τ and a different implementation algorithm, following to the work of Clarke et al. (2013). In this work, it was employed in two separate steps: first, to provide an initial estimate of the glacier thickness, together with the other three models, and second, to calculate a final estimate of the ice thickness with known ice thickness points (from GPR data) constraining the model. In fact, the strength of GlaTE is the integration between model estimates and ground-proof data derived from GPR profiles. GlaTE performs an inversion procedure, constraining the ice thickness results such that they match, with a certain degree of uncertainty, a series of ground-proof data, while following some smoothness requirements, respect the glacier perimeter, and the values at the glacier border are consistent with the terrain slope outside the glacier. The system of equations to be inverted can be summarized into the matrix in Eq. (3):

$$\begin{array}{}\text{(3)}& \left[\begin{array}{c}{\mathit{\lambda }}_{\mathrm{1}}G\\ {\mathit{\lambda }}_{\mathrm{2}}L\\ {\mathit{\lambda }}_{\mathrm{2}}{B}_{\mathrm{gb}}\\ {\mathit{\lambda }}_{\mathrm{3}}{B}_{\mathrm{0}}\\ {\mathit{\lambda }}_{\mathrm{4}}S\end{array}\right]h_{\mathrm{est}}=\left[\begin{array}{c}{\mathit{\lambda }}_{\mathrm{1}}{h}_{\mathrm{GPR}}\\ {\mathit{\lambda }}_{\mathrm{2}}\mathrm{\nabla }{h}_{\mathrm{Clarke}}\\ {\mathit{\lambda }}_{\mathrm{2}}\mathrm{\nabla }{h}_{\mathrm{boundary}}\\ {\mathit{\lambda }}_{\mathrm{3}}\\ {\mathit{\lambda }}_{\mathrm{4}}\end{array}\right]\end{array}$$

where h_GPR is the ground proof GPR-derived ice thickness, h_Clarke is the ice thickness modeled according to Clarke's algorithm, ∇h_bound is the gradient of terrain slope at the boundary of the glacier. The operator G ensures the constraint with GPR data, L with the ice thickness modeled according to Clarke's algorithm, B_gb with the slope outside the glacier, B₀ with the 0 thickness at the boundary, while S is a smoothing constraint. The λ factors are weighting factors and are varied in an iterative manner, in order to give maximum weight to the Clarke model and the smoothness constraint while fitting the GPR data (Grab et al., 2021). The ice density was set at 900 kg m⁻³. The creep factor (or ice softness) A was set to be $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{24}}$ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{3}}$, neglecting its temperature dependence as if the glacier was at 0 °C (Cuffey and Paterson, 2010), since the Rutor Glacier shows temperate conditions (even if this could potentially not be true for the entire glacier) (Cook et al., 2020). The exponent of Glen's flow law was fixed at 3, as considered the best approximation in the absence of data about the ice fabric (Glen and Paren, 1975). The potential for debris presence, which can be added to the model, was not taken into consideration, because both the orthophotos and visual investigation showed no evident thick debris cover in the ablation area of Rutor Glacier. For further details, see Langhammer et al. (2019a), Grab et al. (2021), Schwanghart and Scherler (2014).

The model is open source, it runs in MATLAB environment and is available in an open repository (Maurer, 2022). In the Supplement, the GlaTE model code adapted for this case study is available to ensure full reproducibility. The parameters.txt file in the Matlab folder contains the input parameters used in this study.

2.1.3 OGGM

The OGGM (Open Global Glacier Model) is an open-source collection of algorithms written in Python that provides different insights into glaciers, for example, thickness, meltwater runoff, and future predictions based on climate variations. Its main aim is regional-scale modeling, but the code is modular and can be adapted to work with a single glacier. In this work, the OGGM version 1.6.2 was used, and only the OGGM ice thickness inversion algorithm was employed, which is based on ice flow dynamics and mass conservation (Farinotti et al., 2009; Maussion et al., 2019). The ice flux is computed as Eq. (4):

$$\begin{array}{}\text{(4)}& q=uS=\left(f_{\mathrm{d}}h{\mathit{\tau }}^n+f_{\mathrm{s}}{\displaystyle \frac{{\mathit{\tau }}^n}{h}}\right)S\end{array}$$

where h is the ice thickness, q is the ice volume flux, u is the ice flow velocity, S is the section, which in case of a simplified parabolic section is $=\mathrm{2}/\mathrm{3}\cdot h\cdot \text{width}$, n is the exponent according to Glen's law (=3), τ is the shear stress, f_d is proportional to the ice softness A ($f_{\mathrm{d}}=\mathrm{2}A/n+\mathrm{2}$), f_s is a sliding factor, neglected for simplicity in this run of the model. The flux q in a section is also equal to the mass balance (mass input − output due to precipitation and melting) integrated over the area of the glacier situated above the section considered. This model, in its latest version, works even without the equilibrium assumption, calibrating the mass balance parameters using a geodetic mass balance calibration (Hugonnet et al., 2021). The climate data are taken from the W5E5 dataset (Lange, 2019). During the inversion process, one parameter, the ice softness A, is calibrated against a regional consensus ice volume (Farinotti et al., 2017) and is allowed to vary, in general about one order of magnitude, from the standard value of $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{24}}$ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{3}}$.

Further details and the implementation in the OGGM framework are described in Maussion et al. (2019) and the software is freely available in an open repository (Maussion, 2024). In the Supplement, the OGGM model code adapted for this case study is available to ensure full reproducibility. Among the various folders and files, the .ipnyb Python Jupiter notebook contains a small section about the input and calibrated parameters used in this study.

2.2 Ground-penetrating radar (GPR)

Two different low-frequency antennas were tested to survey the thick ice of the Rutor Glacier. In the 2012 helicopter-based dataset, a GSSI single-frequency antenna with a central frequency of 70 MHz was employed. In the 2022 ground-based survey, a 40 MHz antenna, RIS ONE model with single-frequency antenna configuration, manufactured by IDS, was carried on skis by an operator. The location of the GPR profiles is shown in Fig. 1.

The raw data were processed using the commercial ReflexW software (Sandmeier, 2012), with the following method:

1. Application of a background removal filter (x direction average-trace removal over the entire profiles), to eliminate instrument noise constant in the x direction.

2. Application of a bandpass butterworth filter, from 0 to 150 MHz, to eliminate high frequency noise.

3. Application of a make equidistant-traces filter, to plot 1 trace per meter, especially important for the 2012 survey, since the helicopter was not flying constantly at the same speed.

4. Application of a gain filter called Energy decay: a gain curve in y direction (time) applied on the complete profile based on a mean amplitude decay curve, which is automatically determined (Sandmeier, 2012), to compensate attenuation and geometric spreading in the time-direction.

5. Conversion of the y axis from time to depth, assuming a constant velocity of the electromagnetic signal in the ice of 0.167 m ns⁻¹ (Bohleber et al., 2017). For the helicopter-based survey, we also removed the part of the signal traveling in the air. We did not notice any interference from nearby slopes or features, as in other studies (Church et al., 2018): this is because the helicopter was flying at a low altitude above ground and most of the Rutor Glacier is not surrounded by steep slopes.

6. Manual selection of the ice-bedrock interface with the guidance of model estimations.

The data were not migrated. Some attempts with simplified velocity models were made, without any significant enhancement. Given the very scattered radargrams of this survey, we investigated most of the advanced processing algorithms included in the ReflexW software, in order to find eventual filters improving the radargram clarity. However, simple processing flows, such as the one reported here, gave the most effective results.

Figure 1 The GPR profiles of the Rutor Glacier survey. Star markers indicate the start of the profiles. Numbering from 6 to 10 is based on the numbering of the original helicopter-based dataset (Morra di Cella, 2024). The black glacier margin is that of 2021 used in this study (Macelloni et al., 2022). The pink glacier margin is the 2003 RGI outline (RGI Consortium, 2017).

2.3 2D and 3D visualization of GPR and models together – selection of the bedrock interface

We obtained an overview of the ice thickness by running the three unconstrained models. Then, we needed tools to visualize the models together with the processed GPR data, to perform the selection of the ice-bedrock interface with the model guidance.

The first visualization tool that we employed is RGPR, a software package to analyze and process GPR data (Huber and Hans, 2018). Since it is open-source (it runs in the R environment), we adapted the code to plot the ice thickness estimated by the models together with the GPR profiles. The visual overlap of GPR and models helped us in picking the ice-bedrock interface.

With Paraview software (Ayachit, 2015), we built a 3D reconstruction of surface topography, bedrock topography based on the models, and GPR profiles. The 3D visualization allowed us to directly see where the GPR profiles intersect, helping recognize the ice-bedrock interface as it develops along different profiles. Moreover, visualizing the entire glacier in a 3D environment makes it easier to integrate the user's local knowledge of the glacier in the data interpretation.

Some topographical adjustments were necessary to analyze GPR observations acquired in different years (2012 and 2022). Moreover, the DEM used as input for the models was acquired in 2021. In other words, the models represent the 2021 situation, and the GPR data corresponds to 2012 and 2022.

To visualize together the four ice-thickness models and the GPR sections of the 2012 survey, every model had to be “converted” to 2012, that is, the ice lost from 2012 to 2021 had to be accounted for. In fact, the ice lost in that time frame, especially in the lowest parts of the glacier, was not negligible, as shown in Fig. 2. Since a good 2008 DEM of the glacier surface was available from the regional cartography (Val d'Aosta Region, 2008), the ice loss from 2012 to 2021 was estimated by a simple linear interpolation between the 2008 glacier surface elevation and the 2021 glacier surface elevation, supposing that the glacier melting in those years was constant on average. The calculated ice lost from 2012 to 2021 was added to the model thickness. For comparing the 2022 GPR profile to the models, no correction was performed, since the ice lost from 2021 to 2022 is negligible for the type of comparison performed, and the 2022 profile is situated at a high altitude, where melting is minimal.

To perform these topographical adjustments and analyses, we made sure that all data were reprojected to the same Coordinate Reference System (WGS84, UTM 32 N), using the warp/reproject tool of QGIS software.

Thanks to these visualization tools, the comparison between models and GPR could be done. The ice-bedrock interface was then manually selected (“picked”) in ReflexW, assigning different codes (and colors) to “sure” pickings, “guided by model” pickings, and “probably incorrect” pickings. The “probably incorrect” pickings were selected in those portions of the GPR profiles where it could seem reasonable to “see” an ice-bedrock interface, but considering the broader picture of the glacier, by observing the glacier shape estimated by the models, and based on local knowledge, it was likely a false positive. Those “probably incorrect” pickings served only to show the difficulties in interpreting scattered GPR data, and they were not used in the following analysis of the GlaTE model constrained by GPR data.

2.4 GlaTE model constrained by GPR and sensitivity analysis

After the manual selection of the reflectors in the GPR data, the final GlaTE model, constrained by GPR data, was run. To achieve this, the ice-bedrock interface manually selected on GPR data, based on the 2012 survey, was corrected to 2021, by subtracting the ice lost from 2012 to 2021 from the thickness. The GPR-based ice thickness was given an estimated error of 30 % as GlaTE input parameter, to account for the uncertain nature of the pickings due to the scattered data.

A simple sensitivity analysis was performed on the glaciological parameters A (creep factor [${\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{3}}$]) and gsi (fraction of creep relative to basal sliding [–]). The creep factor A varied from $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{24}}$ to $\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{24}}$, $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{24}}$ being the standard value as suggested by Cuffey and Paterson (2010), while gsi varied from 0.25 to 0.75, with 0.5 being the default value suggested by the model developers. In QGIS, we calculated the average and standard deviation of the various ice thickness models estimated with the different parameters, for both constrained and unconstrained GlaTE models.

3 Results

Figure 2 presents a map of the ice-thickness changes of the Rutor glacier from 2008 to 2021. The ice loss of 100 million m³ does not agree well with the previous GPR-based ice volume estimate of 150 million m³ by Villa et al. (2008), because it would mean that only 50 million m³ still remain in the glacier in 2021, which is improbable, given its size.

Figure 2 The Rutor Glacier (RGI ID: RGI60-11.03039) divided in four categories, according to changes in ice thickness (m) from 2008 to 2021 (“tokyo” color scale, according to Crameri, 2021). The black glacier margin is that of 2021 used in this study (Macelloni et al., 2022). The pink glacier margin is the 2003 RGI outline (RGI Consortium, 2017), and this choice is visible in the yellow areas near the three tongues of the glacier, which did not change in altitude from 2008 to 2021 because they were already free of ice in 2008.

Based on the above observation, we supposed that the previous GPR data by Villa et al. (2008) could have been misinterpreted. In fact, we noticed great difficulties in interpreting our GPR data on the temperate Rutor Glacier, because they were affected by scattering and random noise. Figure 3 shows an example section of the helicopter-based 2012 dataset. An omnipresent clutter zone is situated at 10–50 m of depth, and could be easily misinterpreted for the true bedrock. Fortunately, the true bedrock is visible in the right part of the profile. It submerges below the backscatter zone, suggesting that the former clutter area is probably not a bed return. Elsewhere, possible ice-bedrock interfaces may be spotted and somewhat followed in the GPR section, but the interpretation is far from straightforward.

Figure 3 The GPR section number 7 of the helicopter-based 2012 dataset. The thickness scale was built assuming a constant velocity of the electromagnetic signal in ice of 0.167 m ns⁻¹. The white area is compact firn or ice without water. The ice-bedrock interface is well visible only on the right part of the image.

Therefore, additional information of the bedrock topography was needed, and we retrieved it from ice-thickness models. The four raster maps in Fig. 4 show the model reconstructions of the glacier geometry. The total ice mass estimated by the models was consistently higher than the 150 million m³ of the previous estimate, and it was of about: OGGM = 630, GlaTE = 510, GlabTop2-Py = 580, Original-Glabtop2 = 470 million m³. Average and standard deviation plots are also included in Fig. 4, while a main flowline profile is shown in Fig. 5.

Figure 4 Ice thickness maps of the Rutor Glacier produced with GlaTE, GlabTop2-Py, OGGM, and Original Glabtop2 models without any constrain by ground-proof data, but only with topographic surface data as input.

Figure 5 Main flowline profile. The surface topography is in black, the modelled bedrock is in various colors.

Figure 6 shows an example ice-thickness profile extracted from the models, matching the location of a GPR section. On such profiles, we manually selected the reflections of the ice-bedrock interface. The guidance of the models helped us understand more consciously that the clutter zone identified in Fig. 3 is not the bedrock. The manual selection of the bedrock interface was indicated in Fig. 6 with different colours, to distinguish the “sure” pickings, the probably “incorrect” ones, and those guided by the models. In this way, we selected about 15 % of GPR data as “sure” pickings, 26 % as “guided by models”, and 37 % as “probably incorrect”. Some of GPR data was not picked at all, where the GPR data was too scattered or the models were not consistent between each other. Moreover, a 3D environment was set up in Paraview software to better visualize the GPR profiles as they intersected between each other and the modeled bedrock. This visualization tool enhanced the visual interpretation of the GPR data, and an example is given in Appendix A, in Fig. A7. All the other GPR sections with the model estimations are shown and commented in detail in Appendix A (Figs. A1 to A6). The reader is strongly encouraged to look at those GPR profiles and model estimations, in order to follow the Discussion section more conscious of the subjective interpretation problems we faced.

The bedrock interface selected in the GPR data was then imported into the GlaTE model to perform a constrained run. Figure 7 shows the results of the last GlaTE inversion, performed constraining the model with the manually selected ice-bedrock interface of all the GPR profiles. The total ice volume calculated in this way was about 450 million m³.

Figure 8 shows the results of the sensitivity analysis on the final GlaTE model, in both its unconstrained and constrained version. Maps of average and standard deviation are shown, as well as a profile visualization. By varying the input glaciological parameters A and gsi, the unconstrained version of the model experienced a standard average variation of ±75 millions m³ of ice, whereas the variation of the constrained model was negligible.

Figure 6 An example GPR section with the estimation from the four models overlayed: GlabTop2-Py in red, GlaTE in green, OGGM in blue, Original-GlabTop2 in purple. The manual selection of the ice-bedrock interface, which was possible by comparing the GPR profiles and the models, is displayed in different colours according to their reliability. Purple numbers indicate where this profile crosses other GPR profiles.

Figure 7 The final model of the Rutor Glacier ice thickness obtained by GlaTE model constrained with the GPR data. The pickings of GPR data are shown as differences between GPR ice thickness and GlaTE ice thickness. The final GlaTE model was fairly consistent with the GPR input data.

Figure 8 Sensitivity analysis plots on the unconstrained (top 2 plots) and constrained (middle 2 plots). Profile 7 view of the sensitivity analysis (bottom plot).

In Table 1 we report the ice volume estimates of the different models employed in this study, as well as ice-thickness literature products from Farinotti et al. (2019), Millan et al. (2022) and Villa et al. (2008).

Villa et al. (2008)

Table 1 Ice volume comparison between the different ice-thickness products of this study and literature.

4 Discussion

We propose a methodological approach for jointly interpreting models and GPR data, testing it on the Rutor glacier. We focus on the specific challenge of scattered and unclear GPR data on temperate glaciers and evaluate how ice-thickness models can help interpreting such datasets. First, we compare the four ice-thickness models, then, we discuss their joint interpretation with GPR data, focusing on subjectivity issues.

4.1 Comparison of the four ice-thickness models

The three different models, in their first run without GPR constraints, gave fairly similar estimates of ice volume (OGGM = 630; GlaTE = 510; GlabTop2-Py = 580; Original-GlabTop2 = 470 million m³), with an average of about 540 million m³. In particular, the glacier bedrock shape was similar. Setting the resolution of the models at 20 m avoided the exaggeration of bedrock overdeepenings that resulted from a too-fine resolution. Note that, in the Original-Glabtop2 model, the algorithm automatically resampled the data to 75 m resolution to perform the processing, and the final results was resampled back to 20 m. Observing the distribution of standard deviation, the most uncertain areas are those with overdeepenings. The OGGM estimates were qualitatively similar to those of GlaTE and GlabTop2-Py (which both use the perfect plasticity method), except in some minor parts, despite their differences in the algorithms employed. The OGGM estimated a higher volume, probably due to its mass balance calibration against the consensus estimate (Farinotti et al., 2019), the latter being of about 630 millions m³ (Table 1 and Fig. A9). The consistency between different models, especially in terms of the general shape of the glacier, was already observed in a previous review (Farinotti et al., 2017), which demonstrated that models based on surface topography can provide fairly reliable estimations of ice thickness. They can provide gross estimates of ice volume or future position of lakes in place of glacier overdeepenings (Viani et al., 2016), with less effort in terms of input data and computational time compared to a GPR survey. However, most of the models performed poorly near the margin of the glacier: here, the thickness was generally overestimated compared to the GPR data. We suppose that this fact could be due to the simplification of the glacier geometry and interpolation procedures used by the models, which work better in elongated-shape valley glaciers, but perform worse in wider and more complex glaciers such as the Rutor glacier. The Original-GlabTop2 model was the only one that consistently showed smoother ice-thickness variations near the glacier margins (well visible, for example, in the left part of Fig. A6. The difference between the original and the python implementation of GlabTop2 mainly resides in interpolation and glacier-margin related issues. In particular, the Original-Glabtop2 selects less random points to calculate the ice thickness, and therefore the final model is smoother and the glacier margin cells, which have a fixed thickness of 15 m, weight more in the interpolation process. This observation points out that side aspects of the algorithms, such as the interpolation to the glacier margins, can influence the results in a non negligible way.

4.2 Joint interpretation of GPR and ice-thickness models

The various GPR sections had a different degree of strength of the bedrock return. In some of them, the ice-bedrock interface could be followed relatively easily, such as section 2012 – 8 (Fig. A3). In that case, the four models (GlaTE, GlabTop2-Py, Original-Glabtop2, OGGM) were generally consistent with the ice-bedrock interface depicted by the GPR reflections, demonstrating the overall good quality of the modeled ice-thickness. In most other sections the bedrock reflection was rapidly lost below about 50 m of depth. In general, the signal suffered from scattering, due to the many reflection events distributed throughout the glacier, except in some areas near the margin. No noticeable difference was seen between the two antennas used (40 and 70 MHz). Since the interpretation of these sections was challenging, visualizing models and GPR together was fundamental to retrieving the bedrock topography where it could not be detected by the GPR alone. In particular, in the parts of the profiles where all four models indicated the same ice thickness, we were more confident in selecting that ice thickness as the correct one, provided that it was a reasonable interpretation of the GPR data, although scattered. Mainly, this joint interpretation prevented the mistake of interpreting the first non-reflective layer (white in the GPR sections) as ice and the first reflective zone (scattered black) as bedrock.

A possible source of uncertainty in detecting the correct ice-bedrock interface may also be due to off-nadir reflections, i.e., from valley sidewalls. However, we do not think this issue is critical in this glacier, which is large and has no steep side valley walls.

Specific comments on the interpretation of each GPR section are available in Appendix A. They are different examples of how problematic GPR profiles were interpreted thanks to this methodology, and which mistakes were avoided.

The ice thickness estimated with the second run of the GlaTE model constrained by GPR data is about 450 million m³, which is not too far from the GlaTE estimate without GPR data (510 million m³). Partially, the estimate of the GlaTE constrained model is lower because the models tended to overestimate the thickness near the glacier margin, an issue that was mitigated by the GPR data used in the constrained model. The final estimate of 450 million m³ is almost three times larger than the 150 million m³ previously calculated (Villa et al., 2008). Two literature ice-thickness products were compared to this result, the consensus estimate by Farinotti et al. (2019) and a velocity-based model by Millan et al. (2022). Those models showed consistency between them in terms of volume estimate (both about 630 million m³, even though the bedrock shape was quite different (Fig. A9). The higher volume estimate compared to the GlaTE model is probably due to the fact that these products still rely on the 2003 RGI outline (RGI Consortium, 2017), which is outdated due to the rapid shrinking of the glacier. Those ice-thickness products could be used instead of running the three models, for the same goal of helping to select the correct bedrock interface in the GPR radargrams. It would be simpler than running the three models as performed in this study. However, these regional scale models may have outdated glacier outlines, which are changing rapidly, and therefore thay may provide unreliable estimates in the ablation zone.

4.3 Advantages and limitations of the methodology

Without the help of ice-thickness models, only a small part of the GPR profiles clearly identified the bedrock. This is not an ideal situation to draw a thickness map of the glacier, and it is argued that a similar interpretation problem could have arisen in the 2006 dataset analyzed by Villa et al. (2008), explaining the unexpectedly shallow thickness calculated there. In the GPR profiles, we highlighted in purple the possibly “incorrect” pickings, which mainly represent the clutter zone easily misinterpreted as bedrock, without the model guidance. This clutter phenomenon has been previously reviewed by Rutishauser et al. (2016), which observed that, depending on the glacier, only from 12 to 69 % of the bedrock could be identified. This review was performed on Swiss glaciers, which have climatic conditions similar to the Rutor Glacier. Moreover, better results were achieved only in alpine glaciers where cold conditions are more widespread. Therefore, including modeled ice thickness to guide the GPR interpretation could represent a tool to increase the “pickable”, or “selectable”, regions of scattered GPR sections. This methodology can help especially in the areas where generally GPR is more difficult to interpret: where the ice is thickest, where the bedrock slope is highest, and where there is englacial water (Rutishauser et al., 2016). On the Rutor glacier, the combination of all of these likely challenges the visibility of the bedrock reflection. Another issue of GPR surveys on glaciers that this methodology helped to solve is the limited spatial resolution. The spatial coverage of GPR surveys is limited by survey speed, time and access (e.g. crevasses), leading to discrete, limited sampling of the glacier bed. It is, therefore, possible that the maximum ice thickness remains unknown due to limited survey coverage. Spatial resolution and speed of investigation have always been a compromise. By using the constrained GlaTE model rather than only interpolation, better informed ice-thickness estimates can be made in areas not surveyed by GPR.

The integration of GPR and models seems the most viable option to provide a more reliable estimate of ice thickness than one method alone, especially in the absence of costly boreholes that intercept the bedrock at depth. This is a great achievement, because those models are open-source and require low effort to use, and their reliability and comparability with GPR data has been observed in this and previous research (Farinotti et al., 2017). They also can complement the design of a GPR survey to select the best antenna frequency, based on the expected ice depth, alongside forward modeling to produce synthetic data (MacGregor et al., 2021). Other geophysical surveys, such as Electrical or Seismic tomographies, could also be employed, but to investigate more than 100 m of ice thickness, as needed in this glacier, one would need to set up very long profiles (>500 m) and use powerful seismic or electrical sources. These methods are not as logistically convenient as the GPR.

However, some drawbacks have to be considered.

First, the difficulties in interpreting the GPR data, which represent the most significant source of possible errors and subjectivity; second, inaccuracies in the estimates provided by the different models, which are a simplification of the real system, and rely on estimates of many parameters. However, the most delicate aspect which may limit this methodology is the balance between the reliability of GPR and that of the models. In fact, our approach is a joint interpretation of two different methodologies, in a context where data are of indirect nature and subjective choices must be made. We need the information of the models to better understand the scattered and sparse GPR data, and vice versa, in a sort of manual iterative interpretation process. We stress that, in this difficult context, we do not see the GPR as the main technique, nor the models only as a supporting information, rather, we see that they both support each other. However, the subjective interpretation is potentially prone to human error. In our opinion, the key to minimizing this potential error is to perform simulations on multiple models, display the GPR sections with multiple color scales and in both 2D and 3D environments, and acquire local knowledge about the studied glacier. In this regard, we performed a small, qualitative inter-analyst comparison experiment at a conference in Rome, September 2025 (https://www.gpritalia.it/gpritalia2025/, last access: 9 December 2025). In this experiment, after explaining the same methodology described in this paper to the audience (GPR professionals and academics), we gave the participants a printed version of the GPR profiles, asking them to pick the bedrock of all the profiles, before and after seeing the model estimations. The 10 participants to this activity mostly picked the clutter zone as bedrock (similarly to the clutter zone indicated in Fig. 3 with red arrows). With this type of picking, the total glacier volume would be greatly underestimated, as in the previous work of Villa et al. (2008). Those (unfortunately few) who also reported on their papers the bedrock picking after seeing the models showed that they recognized potential deeper bedrock reflections, closer to the models.

Another source of error is the glacier perimeter, which is an important input of the models, because it was retrieved with a manual observation of the aerial orthophoto. In many glaciers, the perimeter is not straightforward to draw, due to debris hiding the ice (Santin et al., 2023). Moreover, the fact that the DEM and the GPR surveys were carried out at different years required a time interpolation of the ice thickness, to take into account the progressive melting of ice, which introduced another deviation from the real value of ice thickness. The combination of multiple models and methods was key to minimizing the errors; however, the uncertainty remains a main factor in this kind of analysis and it is also difficult to estimate accurately (Aguayo et al., 2024; Shahateet et al., 2025). We observed a standard deviation of about 100 million m³ of ice, which corresponds to an average ice thickness of about ±13 m, by comparing the four unconstrained models; we observed instead a standard deviation of about 75 million m³ by comparing the unconstrained GlaTE model with varied glaciological parameters (sensitivity analysis). These values represent a reasonable estimate of the potential standard deviation in the ice thickness, but they are still biased by the choice of models used and the choice of the parameters to be varied in the sensitivity analysis. Moreover, the definition of standard deviation itself is probably not completely applicable in this case, since it would need a Gaussian error distribution for the ice thicknesses estimated by the various models, which is probably not true. Interestingly, the standard deviation of ice volume calculated during the sensitivity analysis of the constrained GlaTE model is particularly low (about 5 million m³), meaning that the constrained model is quite robust to variations in the input glaciological parameters and tries to adhere to the input GPR data, even if the given GPR data were quite sparse and they were assigned an uncertainty of 30 % in the GlaTE input parameters. It is difficult to say if this very low standard deviation is indicative of real reliability of the GlaTE constrained model; nevertheless, 450 million m³ of ice, distributed as in Fig. 7, represents the best estimate of the Rutor glacier thickness we could obtain with this integrated methodology.

A final and obvious limitation of the study is the absence of any borehole to confirm either the GPR and the model estimated thicknesses. Future research may consider analyzing glaciers where such borehole data is available. Nevertheless, in some decades, due to the ongoing glacier melting, we might directly see the current glacier bedrock in most parts of the glacier. These observations, together with the analysis of DEMs at different years, now possible due to the widespread satellite observations, will provide accurate geodetic mass balance information and surface topography that will deepen the knowledge about ice-flux behaviour in a non-equilibrium state and better calibrate the ice thickness models.

4.4 Future applications and perspectives

For different reasons, both the GPR and ice-thickness models have non-negligible accuracy limitations in reconstructing glacier geometry. The GPR suffers from scattering in temperate glaciers, while models alone have a high uncertainty based on which parameters we choose as input. This study highlights the benefit of combining the two worlds. For GPR applications, the use of models can help in designing the survey and selecting the proper instrumentation (e.g., antenna frequency) based on the estimated geometry of the glacier. Models can also provide a picture of the glacier in those zones inaccessible to GPR for logistics limitations or time constraints. For modeling studies, including GPR surveys in their workflow can improve their estimations, especially for those glaciers with complex and difficult-to-model geometries. Thus, applying the workflow presented here to other glaciers with GPR data available could help calibrate regional models more accurately, or at least raise awareness of the potential uncertainties of the models. The Supplement can be used to reproduce the same workflow on other glaciers.

Future integration of GPR data in glacier models could include not only ice thickness, but other features of the glacier which are detected by GPR and are very difficult to predict with models, such as the distribution of cavities and the distinction between cold and warm ice zones inside a glacier. The latter, according to Reinardy et al. (2019) and Comiti et al. (2019), plays an important role in regulating the sediment transport at the glacier base and can be retrieved from the distribution of reflections (due to englacial water) in the GPR sections. This could enhance the link between glacier geometry and hydrological output in a glacier model. The clutter itself in the radargram can be used to infer glacier features (Scanlan et al., 2020).

From a local perspective, a more robust reconstruction of the Rutor glacier compared to previous estimates will be helpful for ongoing studies on the water mass balance and sediment transport within the Rutor basin. Since the Rutor basin hosts many proglacial lakes, the bedrock topography map reveals the possible position of future lakes, in place of the overdeepenings of the bedrock topography (Fig. A8). Their location is consistent with the estimates of a previous study (Viani et al., 2020). The shrinking of Rutor Glacier is speculated to occur mostly in the following decades, given its volume of about 450 million m³ and its loss of 100 million m³ from just 2008 to 2021. After a few decades, little ice is expected to be still stored in the Rutor Glacier.

5 Conclusions

Investigating glacier substructures with GPR may be challenging in temperate glaciers, where widespread water content and debris cause signal scattering, making it difficult to distinguish the ice-bedrock interface. The Rutor Glacier had already been surveyed with GPR in the past but, due to these interpretation difficulties, was estimated to have a very small ice thickness, of about 17.5 m on average. This estimate proved to be wrong after observing that, in the 2010–2020 decade, the Rutor Glacier lost more than 20 vertical meters in $\mathrm{1}/\mathrm{4}$ of its area, while reducing its area only by a fraction.

The analysis of two new GPR datasets from 2012 (helicopter-based) and 2022 (ground-based) confirmed the difficulty of reliably detecting the ice-bedrock interface. Therefore, the Original-Glabtop2 and the open-source GlabTop2-Py, GlaTE, and OGGM models were tested, to understand how they could support the interpretation of difficult datasets acquired on temperate glaciers. First, these models were run with only the glacier surface topography as input. Then, their estimated thickness was overlapped with the GPR sections, providing substantial help in manually selecting the ice-bedrock interface. In particular, this methodology avoided the misinterpretation of englacial water-rich areas as the ice-bedrock interface. Finally, a second run of GlaTE produced an ice-thickness model of the Rutor Glacier constrained by GPR data.

Analyzing the ice-thickness variations among models and by varying the glaciological parameters of the GlaTE model, we recognized that the uncertainty in both GPR and models is a major factor and it is difficult to measure a true value of ice thickness in a temperate glacier.

A preliminary run of two or three ice-thickness models, such as the ones tested in this study, is advised before carrying out a GPR survey on a glacier. The ice-thickness models, in combination with the GPR, are the most cost-effective way to represent, with a certain degree of uncertainty, the glacier bedrock geometry. This is in line with the philosophy behind the GlaTE model, built to constrain a topographical-data-based algorithm with “ground-proof” GPR data. An effort was made to provide supplementary materials which can be used to reproduce the same workflow on other glaciers. Local studies can benefit from more accurate glacier geometry estimates, as well as regional studies, which can be calibrated more effectively.

Appendix A: GPR Sections

This appendix is devoted to all the GPR sections analyzed in this work, to show the interested reader what the data looked like, which were the difficulties, and in which cases the three models' predictions were useful in avoiding clear misinterpretations. In any case, the interpretation subjectivity is high, and analyzing the openly available GPR dataset with specialistic software is advised. The first 5 figures represent the 5 GPR sections of the 2012 helicopter-based survey; the last figure is the merging of all (subsequent) sections of the 2022 ground-based survey. Purple numbers on top of the figures indicate where that profile crosses other GPR profiles. See Fig. 7 for the location of the GPR sections.

Additional figures at the end of the appendix show: a 3D visualization of models and GPR in Paraview (Fig. A7); a map of the bedrock and surface topography as in the constrained GlaTE model, visualized in contour lines (Fig. A8); a comparison with literature ice-thickness products (Fig. A9).

Figure A1 GPR section 2012 – 6. This section is situated at the top of the glacier, and in its central part (500–1000 m of Distance coordinate) the three models estimate an overdeepening of over 150 m. The GPR reflections are reasonably clear until 50–70 m of depth. In the right part of the picture, the models follow closely the GPR reflections; however, they are quite imprecise just near the margin, where they do not draw correctly the bedrock shape where the ice is very thin. Based on the estimates of the models, the analyst did not trust the sparse reflections at 25–30 m of depth as an ice-bedrock interface but acknowledged that, in the right part, the deepening reflections continue to go deeper leftward to around the center of the image.

Figure A2 GPR section 2012 – 7. This section was already shown and discussed in Fig. 6, but it is possible to add a consideration on the steep overdeepening “seen” by the GlabTop2-py model. The creation of deep overdeepenings was recognized to be an artifact of such models, especially when using a too-fine input DEM [Maurer and Grab, personal communication]. Also, the other models estimate a high ice thickness at that point, but they seem more realistic. Unfortunately, the GPR was of little help in that region.

Figure A3 GPR section 2012 – 8. This GPR section was probably the most clear of the dataset. The ice-bedrock reflector can be followed to great depths, although the clarity is far from ideal. This section proves that the models can reasonably estimate the ice thickness of the Rutor Glacier, because they match with GPR data where the latter are clear. In the left part of the image, which corresponds to the Eastern tongue, the models are not shown because they were cut with the 2021 margin. From 2012 to 2021, the glacier experienced a notable retreat in the Eastern tongue.

Figure A4 GPR section 2012 – 9. This GPR section was probably the least clear. Its path runs along the elongated overdeepening of the top of the glacier (Fig. 7). In the left part of the picture, the reflections deepen very fast, and this is displayed also by the models, but after that, it was impossible to retrieve any other reflections attributable to the ice-bedrock interface. The role of the models, in this case, was very important, to avoid thinking that the darker reflections seen in the center-right part of the image at 30-50 m depth are the bedrock interface, while probably it is due to the higher englacial water content, or debris (Forte et al., 2021).

Figure A5 GPR section 2012 – 10. This GPR section was generally clear in the right part of the image, and also showed consistency between GPR and models. All the profile has an interface around 30–40 m interpreted as englacial water or debris (Forte et al., 2021). The left part of the picture was more problematic, because the GPR, although not clear, seems to suggest a deeper bedrock interface than the models. Probably, this is due to a peculiar location of this survey line, possibly running along a very high-sloping bedrock (Fig. 7). Such high-slope bedrock areas are known to disrupt GPR measurements and they could be common near the overdeepenings of Rutor Glacier since they are also evidenced in many locations in the proglacial zone (which was formerly occupied by ice during the past ice ages).

Figure A6 GPR section 2022. This GPR section was acquired with an antenna (40 MHz) different from that used for the previous sections shown. It also was ground-based and not helicopter-based. Notwithstanding the lower frequency and being ground-based, the GPR reflections were not particularly clearer, showing a physical limit of the technology in the presence of widespread englacial water or debris (Forte et al., 2021). The models did not perform well near the margin, with the exception of Original-Glabtop2. as observed also in the other GPR sections. However, probably the models offer a reasonable estimate of the bedrock interface in all the left-central part of the image. The right part seems more problematic because the GPR signal is completely lost, however, its greater depth is supported by the GPR signal texture: compare e.g. the image at 100 m depth at Distance =1000 and 2500 m. At distance =2500 m, the image shows a granular texture similar to where there is ice and englacial water, while at Distance =1000 m the texture is transparent. Similar considerations were used to interpret the GPR sections by Forte et al. (2021).

Figure A7 3D visualization in Paraview of GPR profiles and bedrock as estimated by the three models. This 3D visualization, even if difficult to render in a 2D image, was helpful to see the GPR profiles at their intersections, where their bedrock interface must match.

Figure A8 Bedrock and surface contour lines as in the constrained GlaTE model.

Figure A9 Comparison between the ice-thickness product of the constrained GlaTE model and two literature products: an ice-thickness map by Millan et al. (2022) and the consensus estimate (Farinotti et al., 2019).