1. Introduction
Estimating glacier volumes is crucial for assessing climate change effects, as glaciers are key indicators of climate change and directly influence sea level, water availability and ecosystems (e.g. Meier, Reference Meier1984; Dowdesewell and others, Reference Dowdeswell, Hagen, Björnsson, Glazovsky, Harrison, Holmlund, Jania, Koerner, Lefauconnier, Ommanney and Thomas1997). Direct measurements of glacier thickness using Ground-Penetrating Radar (GPR) are limited by the fact that they are time-consuming, while boreholes are expensive and provide only point measurements. In addition, these investigations are not practical for hundreds or thousands of glaciers. It is therefore necessary to use various models and algorithms that combine physical, empirical, and remote sensing approaches. These models range from basic empirical relations to sophisticated numerical simulations, each with their unique strengths and limitations depending on the glacier type and data availability.
Mass-conservation approaches estimate glacier thickness by balancing the input from snow accumulation with the output due to ice flow and melting. Several modelling algorithms have been implemented on the principles of mass balance and ice flux dynamics. Glacier thickness can be estimated through ice-flow dynamics, where the thickness is derived from the relationship between ice velocity, surface slope and basal shear stress. In this approach, Glen’s flow law (Glen, Reference Glen1955), which describes ice deformation under pressure, plays the central role. This law, extended by Nye (Reference Nye1970), relates ice strain to the applied stress, allowing for the calculation of thickness when ice velocity and slope are known. One of the most widely used mass-conservation models was developed by Farinotti and others (Reference Farinotti, Huss, Bauder, Funk and Truffer2009), which integrates surface mass balance and ice flux divergence to calculate thickness distribution. This method requires surface velocity data often derived from remote sensing, and digital elevation models (DEMs) to infer the ice thickness along specific glacier profiles. Improvements to this model include the Bayesian Ice Thickness Estimation (BITE) model (Werder and others, Reference Werder, Huss, Paul, Dehecq and Farinotti2020), which uses Bayesian inference to combine prior knowledge with observational data, providing uncertainty estimates along with the thickness predictions. Brinkerhoff and others (Reference Brinkerhoff, Aschwanden and Truffer2016) used a similar approach to make statistical estimates on glacier topography. Another ice thickness modelling approach relies on the shear-stress calculation (Haeberli and Hölzle, Reference Haeberli and Hölzle1995), in which the mass balance is not considered (Farinotti and others, Reference Farinotti2017).
Numerical inversion techniques have become increasingly popular for estimating glacier thickness by using surface velocity data to back-calculate the ice thickness distribution. This is often done by minimizing the discrepancy between observed and modelled values, as described in Morlighem and others (Reference Morlighem, Rignot, Seroussi, Larour and Dhia2011). Another widely used method is volume-area scaling, which provides a simplified way to estimate glacier volume based on surface area measurements. This method is particularly useful for large-scale assessments, where direct thickness measurements are unavailable. The principle behind volume-area scaling is that glacier volume can be approximated as a power-law function of the glacier’s surface area. Bahr and others (Reference Bahr, Pfeffer and Kaser2015) extensively reviewed this technique, highlighting its effectiveness in estimating glacier volume at regional and global scales. However, all these methods have limitations, particularly when applied to debris-covered or highly irregular complex glacier systems, where scaling relationships may introduce biases (Banerjee and others, Reference Banerjee, Patil and Jadhav2020).
The advent of remote sensing technologies has significantly advanced glacier thickness estimation. In fact, satellite-derived and air-borne based data, such as DEMs and ice velocity from synthetic aperture radar (SAR) and optical imagery (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022; Piermattei and others, Reference Piermattei2024) have become critical inputs for many thickness and ice velocity models. These datasets enable large-scale and high-resolution glacier monitoring, allowing thickness models to be applied over large regions. Unfortunately, while repeated airborne and satellite-based altimetry provide accurate surface elevation change data, which can be directly used to infer ice thickness variations over time, estimating the actual ice thickness and volume of frozen bodies at different scales remains a challenge. An additional issue is the presence of different and often coexisting frozen materials, ranging from fresh snow to firn and different types of ice, all with different densities and rheology. None of these materials are pure but are often mixed with debris (Santin and others, Reference Santin, Roncoroni, Forte, Gutgesell and Pipan2024), may contain voids of various sizes, and water (Paterson, Reference Paterson1994).
Glaciers also have different thermal regimes depending on their geographical location, altitude and thickness. According to their thermal regime, glaciers are generally divided into cold, warm and polythermal. Specifically, polythermal glaciers are characterised by the coexistence of both warm (at 0°C) and cold (below 0°C) ice. Typically, it is the thicker, higher ice in the accumulation area that is warm, and it is the snout, lateral margins and surface layer of the glacier that are below the pressure melting point (Glasser, Reference Glasser, Singh, Singh and Haritashya2011).
Distinguishing between cold and warm ice at different scales is critical to understanding glacier dynamics, as the thermal properties of ice influence glacier flow, meltwater production and response to climate evolution. Cold ice exists below the pressure melting point, while warm ice is at or near the melting point and contains liquid water within its structure. It is debated if even cold ice can contain liquid water (Ball, Reference Ball2009). Several techniques are used to distinguish between cold and warm ice, ranging from field-based measurements to large-scale remote sensing technologies.
At the local scale, thermistor probes inserted into boreholes provide direct temperature measurements of glacier ice, allowing precise mapping of thermal profiles (Paterson, Reference Paterson1994; Cuffey and Paterson, Reference Cuffey and Paterson2010), but with the inherent and insurmountable limitation of providing only local measurements.
At the regional scale, remote sensing techniques can be exploited. SAR is widely used for large-scale monitoring of glaciers and ice sheets because warm ice, being softer and more prone to deformation, exhibits different velocity patterns compared to cold ice. SAR can capture surface velocity changes, which indirectly help in identifying thermal properties of the ice (Rignot and others, Reference Rignot, Jezek and Sohn1995; Joughin and others, Reference Joughin, Kwok and Fahnestock1996; Tedesco, Reference Tedesco2015). Thermal differences between cold and warm ice can also be inferred using satellite altimetry data combined with optical imagery because surface elevation changes due to melt or internal deformation provide clues about the thermal state of the glacier (Flament and Rémy, Reference Flament and Rémy2012). The main limitation is that the information is limited to the glacier surface, and the transition between cold and warm ice, as well as the contact between ice and bedrock, cannot be investigated.
At the intermediate scale, geophysical techniques, and specifically GPR is frequently used to detect the boundary between cold and warm ice. Cold ice has lower dielectric permittivity compared to warm, water-rich ice and such liquid water produces diffuse scattering that can be detected and mapped on GPR data (Murray and others, Reference Murray2000; Pettersson and others, Reference Pettersson, Jansson and Blatter2004). Moreover, quantitative estimates can be obtained from GPR, inferring the water content (Gusmeroli and others, Reference Gusmeroli2010), the density variations (Forte and others, Reference Forte, Dossi, Colucci and Pipan2013) and the dielectric properties of the ice (Bradford and Harper, Reference Bradford, Harper and Brown2009; Liu and others, Reference Liu, Takahashi and Sato2014). The main problems are related to GPR data interpretation as the electromagnetic (EM) facies of cold and warm ice are not always so clear (Forte and others, Reference Forte2021), but dedicated signal attribute analyses can, at least partially, overcome such problems (Gutgesell and Forte, Reference Gutgesell and Forte2024a).
Ice thermal modelling algorithms have been implemented in recent years and applied at various scales and for different purposes. Such algorithms are essential for understanding the thermal dynamics of ice sheets, glaciers and sea ice and allow to simulate the temperature distribution within the ice and its interaction with underlying water, atmosphere, and bedrock. There are different approaches including Finite Difference (e.g. McCarthy and others, Reference McCarthy and Waddington2010), Finite Elements (e.g. Franca and others, Reference Franca, Hauke and Masud2006) or coupled thermal and flow dynamics (e.g. Hindmarsh and Le Meur, Reference Hindmarsh and Le Meur2010). Other models are based on empirical or semi-empirical equations (Clarke, Reference Clarke2005) or coupled climate-ice models which integrate different processes with ice thermal dynamics, providing a holistic view of ice behaviour in response to climate change (e.g. Bitz and others, Reference Bitz, Holland, Weaver and Eby2001). Each of these algorithms and methods has its strengths and weaknesses, and the choice of which to use often depends on the specific characteristics of the ice being modelled and the available data and, even more importantly, the desired scale of the modelling. A comprehensive up-to-date review of thermal modelling of ice is beyond the scope of this paper and can be found in Aljuneidi and others (Reference Aljuneidi, Heine, Rodriguez and Boetcher2024).
While there are several scientific papers comparing ice thickness estimates obtained from GPR and different modelling approaches (e.g. Farinotti and others, Reference Farinotti2017; Vergnano and others, Reference Vergnano, Franco and Godio2024), there are only a few papers investigating the thermal structure of polythermal glaciers using integrated modelling and GPR techniques (Delcourt and others, Reference Delcourt, Van Liefferinge, Nolan and Pattyn2013; Wilson and others, Reference Wilson, Flowers and Mingo2013). As far as we know, there are no studies that combine all these elements. Therefore, the motivation of this work is to address the following scientific objectives: (i) To make a rigorous comparison between ice thickness inferred from two GPR datasets and estimates derived from different modelling approaches and empirical laws. (ii) To compare warm and cold ice areas detected by GPR datasets with 2D thermal modelling conducted along the same path of some selected GPR profiles.
All analyses and modelled data are based on two nearby polythermal local glaciers in Southwest Greenland, as detailed in the following section, but we believe that the results are of general relevance and can be useful also in different geographical and climatic conditions.
2. Study area
The two glaciers investigated in this study are among a group of 100 glaciers (Rastner and others, Reference Rastner, Bolch, Mölg, Machguth and Paul2012) located about 30 km northeast of the town of Sisimiut (66°55’ N, 53°40’ W) in Southwest Greenland (Qeqqata Municipality) (Figure 1a). Here and in a previous work (Marcer and others, Reference Marcer2017), these are referred to as the Sisimiut Glaciers (Figure 1b) and cover an area of approximately 85 km2. The Sisimiut Glaciers are located more than 100 km from the Greenland Ice Sheet (GrIS) and are far from other glaciated areas of coastal West Greenland, with the closest being Maniitsoq (Sukkertoppen) 60 km to the south and Qeqertarsuaq (Disko Island), 250 km to the north. Similarly to other glaciers outside the GrIS in West Greenland, they showed a substantial decrease in area (−20 %) and surface elevation (−0.43 m a−1) during 1985–2020, while their equilibrium line altitude shifted 97 m upwards in the same period (Securo and others, Reference Securo2024).
Figure 1. Greenland (Kalaallit Nunaat) (a), Sisimiut glaciers in West Greenland (b) and mt. Aqqutikitsoq with the recently (30 Aug. 2023) installed AWS and the GPR data analysed in this study: western (in blue, 2014) and eastern (in red, 2024) Aqqutikitsoq glaciers (WG and EG, respectively) and GPR profiles (c). Elevation data are retrieved from arcticdem mosaic (Porter and others, Reference Porter2023), glaciers divides are retrieved from Randolph glacier inventory (RGI Consortium, 2023) and water–land vector masks are retrieved from Greenland (Moon and others, Reference Moon, Fisher, Stafford and Thurber2023).
None of the Sisimiut Glaciers has an official name, but one of them has already been studied in Marcer and others (Reference Marcer2017) and is currently monitored with geodetic measurements on an annual basis, since 2021 (unpublished data). Several of these glaciers, including the two investigated in this paper, are frequently crossed by local people in winter, both on skis and snowmobiles, and rarely in summer. Part of the Sisimiut Glaciers meltwater runoff ends into a lake (Isortuarsuup Tasersua) that is used by the local hydropower plant.
In the southern area of Sisimiut Glaciers, south of Isortuarsuup Tasersua and about 1 km from the end of Kangerluarsuq Ungalleq (second fjord), is Mount Aqqutikitsoq (1448 m a.s.l.). The massif is surrounded by several glaciers on its north, south and east sides. Two of these are the focus of this study, named here according to their relative position throughout the text: Aqqutikitsoq Western Glacier (WG) and Eastern Glacier (EG). Both are among the 10 largest glaciers in the area.
WG (RGI v7.0, Region 05, ID 03155) is a valley glacier of 3.3 km2 with a median elevation of ∼950 m a.s.l. (RGI Consortium, 2023). The glacier has two tongues terminating to the east, one of which ends into a small proglacial lake. A third glacier tongue further west is part of the same glacier complex but is unnamed, not included in this study (Figure 1c), and flows in the opposite direction. Marcer and others (Reference Marcer2017) calculated a surface elevation change of −0.60 ± 0.11 m a–1 for the WG from 1985 to 2014.
EG (RGI v7.0, Region 05, ID 03158 and 03164) is part of a larger glacier complex of 11.8 km2 that flows in different valleys, all within the same catchment that ends up into the Isortuarsuup Tasersua. EG area is 7.2 km2 with a median elevation of ∼1000 m a.s.l. (RGI Consortium, 2023). The surrounding topography is less steep than its western counterpart and consists to a lesser extent of avalanche prone slopes above 30 degrees. WG and EG differ in aspect and flow direction, toward West and East, respectively. However, both are of comparable size and can be considered typical among Sisimiut Glaciers larger than 1 km2.
3. Methods
This section describes the methodological details of GPR data acquisition and processing (Section 3.1) and signal attribute analysis (Section 3.2). We also give details about the inference made on the temperature data of the study area (Section 3.3), as this is one of the essential inputs to both the thickness/volume and thermal modelling, which are detailed in Sections 3.4 and 3.5, respectively.
3.1. GPR
GPR is a non-invasive near-surface geophysical technique based on the propagation of high frequency EM waves (typically in the 10 MHz—2 GHz range) in the subsurface. Since signal propagation is mainly affected by the EM properties of the soil, GPR can be used to investigate subsurface structures and materials with varying degrees of resolution and penetration (Jol, Reference Jol2009) in a wide variety of applications in various fields, including archaeology, geology, engineering, hydrology, glaciology and many others (Daniels, Reference Daniels2004). This is mainly due to the GPR’s greater versatility, resolution and acquisition rates as compared to other geophysical techniques. Most GPR surveys are carried out using ground-coupled systems, where the antennas are placed either directly on the ground surface or just a few centimetres above it, and are moved across the survey area by hand or vehicles. However, airborne GPR surveys are preferable in logistically challenging and potentially dangerous areas, as well as for surveying large areas (Rutishauser and others, Reference Rutishauser, Maurer and Bauder2016), potentially covering several tens of kilometres per hour. Similar data collection rates can be achieved for ground-based surveys when the instrument is mounted on vehicles such as cars (Liu and others, Reference Liu2021), all-terrain vehicles (Novo and others, Reference Novo, Lorenzo, Rial and Solla2012) or snowmobiles (Holbrook and others, Reference Holbrook, Miller and Provart2016). GPR is an excellent geophysical method for glaciological investigations in both glacial and periglacial environments for a wide range of applications (e.g. Arcone and others, Reference Arcone, Lawson and Delaney1995; Bradford and others, Reference Bradford, Harper and Brown2009; Forte and others, Reference Forte, Basso Bondini, Bortoletto, Dossi and Colucci2019) because it exploits the generally low electrical conductivity of frozen materials to reach penetration depths that are not possible for most geologic materials.
We use two distinct GPR data sets, namely the ‘western’ and the ‘eastern’ ones, collected on EG and WG, respectively. The former dates from 2014 (August) and consists of 19.3 km of GPR data collected with a 50 MHz Malå unshielded rough terrain antenna carried by hand and connected to a ProEx GPR system. The latter was collected in 2024 (February) and consists of 25.3 km of GPR data acquired with a 100 MHz Malå shielded antenna carried by a snowmobile. Details of both surveys are summarised in Table 1. More information about the 2014 survey can be found in Marcer and others (Reference Marcer2017). Even though the two surveys are very close to each other (Figure 1), there are relevant differences in the antenna type and frequency used, the year, the season, the logistics of acquisition, as well as the spatial density, making the two datasets an excellent test for the comparison of GPR versus ice thickness and thermal modelling algorithms results.
Table 1. Synthesis of GPR surveys and interpretation details
GPR surveys of Greenland glaciers outside the GrIS are not common, although a few examples surveyed during the last decades can be found for the Arcturus and Schuchert Glaciers (E-Greenland; Citterio and others, Reference Citterio, Mottram, Hillerup Larsen and Ahlstrøm2009), the Mittivakkat Glacier (SE-Greenland, Yde and others, Reference Yde2014), the Qasinnguit Glacier (SE-Greenland, Abermann and others, Reference Abermann, Van As, Petersen and Nauta2014), the Lyngmarksbræen Ice Cap (W-Greenland, Gillespie and others, Reference Gillespie, Yde, Andresen and Citterio2023) and the Qaanaaq Ice Cap (NW-Greenland, Lamsters and others, Reference Lamsters2024).
Both the WG and EG GPR datasets were processed with the same processing flow, but setting different parameters considering the different equipment and setting of the two surveys. In particular, we applied a processing flow including: drift removal (zero-time correction), bandpass filtering (corner frequencies of the filter equal to 5–20–100–200 MHz and to 15–30–200–250 MHz for the 2014 and 2024 surveys, respectively), exponential amplitude recovery, topographic (static) correction and f–k time migration. The EM velocity used for both depth conversion and migration was set equal to 0.168 m ns–1 according to Marcer and others (Reference Marcer2017). This value is typical for pure cold ice (e.g. Pettersson and others, Reference Pettersson, Jansson and Holmlund2003) and the assumption of homogeneous velocity is similar to that made in several other glaciological studies (Saetrang and Wold, Reference Saetrang and Wold1986; Melvold and Schuler, Reference Melvold, Schuler, Hauck and Kneisel2008; Yde and others, Reference Yde2014). Considering that in polythermal glaciers cold and warm ice coexist, in the latter the water present between the ice crystals reduces the EM propagation velocity (Bradford and Harper, Reference Bradford and Harper2005). However, estimates of EM velocities in polythermal glaciers suggest that the differences between cold and warm ice are relatively small. For instance, Bradford and others (Reference Bradford, Harper and Brown2009) analysing a multi-offset GPR data set on a glacier in Alaska, concluded that the mean velocity for the cold ice was 0.170 m ns–1, while for the warm ice was 0.160 m ns–1, (i.e. a difference of about 6%) both with large local variations and with uncertainties of about 2% when using reflection tomography algorithms on multi-fold data and between 5–10% when using migration velocity analysis on single-fold (i.e. common-offset) data, as in the present and most of the cases. Similar results are reported by Murray and others (Reference Murray, Booth and Rippin2007) for multi-fold GPR data on two glaciers: one in the European Alps and the other in the Svalbard Islands. They found mean velocities of 0.1624 (± 0.001) and 0.1701 (± 0.007) m ns–1 for the drier shallower ice layer, and of 0.1506 (± 0.02) and 0.1619 (± 0.008) m ns–1 for the deeper layer containing some liquid water, for the two glaciers, respectively. Therefore, the error in the GPR depth conversion and ice thickness estimates introduced by keeping the EM velocity constant and equal to 0.168 m ns−1 is of the same order of magnitude as the uncertainties in the EM velocity analysis techniques, especially for common-offset GPR data sets and does not exceed a few percent. We did not consider the snow cover as its thickness was zero during the acquisition of 2014 data set, and was negligible (the mean value in the accumulation zone was 1.8 m) with respect to the ice thickness in the 2024 dataset.
3.2. GPR attributes
Signal attributes have been originally developed in geophysics to improve the interpretation of reflection seismic data. An attribute is basically ‘any measurement derived from (seismic) data’ (Sheriff, Reference Sheriff2002). Such a broad definition encompasses an incredible number of attributes, based on a variety of algorithms and computational strategies that play an extremely important role in the interpretation and analysis of seismic data (Chopra and Marfurt, Reference Chopra and Marfurt2005; Fomel, Reference Fomel2007). More recently, various signal attributes have been used in GPR data analysis, interpretation and quantitative data extraction, even extending their meaning and computational strategies (Roncoroni and others, Reference Roncoroni, Forte and Pipan2024). They have demonstrated their effectiveness for several glaciological applications (Forte and others, Reference Forte, Dalle Fratte, Azzaro and Guglielmin2016; Zhao and others, Reference Zhao, Forte, Colucci and Pipan2016; Church and others, Reference Church, Bauder, Grab and Maurer2021) since they are able to improve the interpretation of various glaciological structures.
Here we focused on amplitude-, frequency-, and phase-related attributes to better image and understand the boundaries and the geometries of the GPR horizons. In addition, because each geological material has unique physical, chemical and rheological properties, they create a distinctive pattern of reflections, diffraction, attenuation and spectral characteristics, often referred to as EM facies (Gutgesell and Forte, Reference Gutgesell and Forte2024b), similarly to geological facies (Moore, Reference Moore, Longwell, Moore, McKee, Müller, Spieker, Wood, Sloss and Krumbein1949). EM facies thus represent the cumulative appearance (or signature) of a material on a GPR section in response to the propagation of the EM wave.
GPR attribute calculation, interpretation and facies assessment were conducted using Petrel software (Schlumberger), which was originally designed for reflection seismic but can be easily adapted to GPR. Further details on the calculation, meaning and use of GPR attributes in glaciology are beyond the scope of this paper and can be found in Zhao and others (Reference Zhao, Forte, Colucci and Pipan2016).
3.3. Temperature
We used the Copernicus Arctic Regional Reanalysis (CARRA) (Schyberg and others, Reference Schyberg2020) surface meteorological variables forecast to calculate the 1991–2020 mean annual air temperature (MAAT) on both glaciers. We averaged CARRA 24 h daily values into annual averages for the glacier areas using zonal statistics and RGI polygons (RGI, 2023). The 1991–2020 MAAT for Aqqutikitsoq Glaciers was −7.8 ± 1.4°C, ranging from −10.5°C in 1992 to −3.8°C in 2010. In the first 10 years of the reference period, MAAT was −8.5°C, in the second −7.0°C and in the third −7.2°C. The mean summer air temperature (MSAT) and winter air temperature (MWAT) are 2.4 ± 1.1°C and −16.9 ± 2.6°C, respectively. February is the coldest month (−17.9°C) and July is the warmest (3.9°C).
An automatic weather station (AWS; Figure 1c), was installed at the end of August 2023, and is now available at a suitable location for the studies of the Aqqutikitsoq Glaciers. The AWS is located at an elevation of 840 m a.s.l. and a linear distance of 3 and 7 km from the western and eastern glaciers, respectively. We used air temperature from the AWS from 30 August 2023 to 29 August 2024 to assess the accuracy of CARRA data at the AWS location. The difference of MAAT in the same period and location between the measured data and the reanalysis product is −0.01 ± 2.45°C, although the CARRA cell orography value in the AWS location is lower than at the station, at 671 m a.s.l.
3.4. Thickness/volume modelling
To obtain information on ice thickness distributions and therefore ice volumes of glaciers independently of GPR data, we used: modelling, data already calculated from dedicated repositories, and global-scale power–law relations. To apply that relations the areas of the Aqqutikitsoq East and West Glaciers are manually digitised based on available Google Earth imagery, acquired via Pléiades satellites (CNES Airbus) in the area in July 2019. The glacier surface topography is retrieved from ArcticDEM Mosaic (Porter and others, Reference Porter2023). In the study area, this stacked DEM with 2 m resolution is based on ArcticDEM Strips (Porter and others, Reference Porter2023).
As far as ice thickness modelling we used the GlabTop—Glacier bed Topography—(Linsbauer and others, Reference Linsbauer, Paul and Haeberli2012) algorithm. This shear-stress-based model integrates various factors, such as topographic information, geometrical data, and climatic conditions. We used GlabTop2 here (Frey and others, Reference Frey2014), which builds on the foundation of GlabTop and offers improvements in usability. The main improvement of GlabTop2 is that it avoids the tedious process of manually drawing branch lines.
GlapTop2 for the ice thickness (h) uses the following equation (Frey and others, Reference Frey2014):
\begin{equation}h = \tau /\left( {f\rho gsin\left( \alpha \right)} \right),\end{equation}where τ is the average basal shear stress along the central flowline of the glacier [kPa] (Haeberli and Hölzle, Reference Haeberli and Hölzle1995), g is the gravitational acceleration [9.81 m s−2], α is the surface slope [°], ρ is the density of the ice [900 kg m−3], and f is a dimensionless shape factor to be constant at 0.8 as a default value for valley glaciers. To evaluate the sensitivity to f we specifically calculated it by averaging the values separately obtained for EG and WG using the following equation:
\begin{equation}{f_i} = \frac{2}{\pi }arctan\left( {\frac{{{w_i}}}{{2{h_i}}}} \right)\!,\end{equation} where 
${w_i}$ is the width of the i glacier cross-section (using the glacier margins from RGI Consortium, 2023) and 
${h_i}$ is the average ice thickness along it, as estimated by GPR data. We set 18 and 21 cross-sections for the EG and WG, respectively (Table 1S).
In addition to the above-described modelling, we download the thickness of the two study glaciers directly from Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022), which provides data for all the glaciers in the RGI repository (RGI Consortium, 2023).
To further extend the comparison of ice thickness and volume of the glaciers, we selected some global-scale power–law relationships. We applied the relation originally proposed by Macheret and others (Reference Macheret, Cherkasov and Bobrova1988) and validated by Bahr (Reference Bahr1997), with the two crucial parameters (i.e. the power-law coefficient, c and the exponent, γ) set according to different authors: Chen and Ohmura (Reference Chen and Ohmura1990); Radić and Hock (Reference Radić and Hock2010); Adhikari and Marshall (Reference Adhikari and Marshall2012); Grinsted (Reference Grinsted2013). These two parameters are used in equation (3) to estimate the ice volume (V) from the glacier area (A).
\begin{equation}V = c{A^\gamma },\end{equation}3.5. Thermal modelling
In order to obtain an estimate of the glaciers’ thermal regime, independent from the GPR data, we used PoLIM (Polythermal Land Ice Model) a 2D higher-order thermo-mechanical flow band model (Wang and others, Reference Wang, Zhang, Xiao, Ren and Wang2020). PoLIM is a model designed to simulate the behaviour and dynamics of polythermal glaciers and ice sheets. It accounts for the thermal stratification of ice, allowing it to simulate both the cold ice that is below freezing and the temperate ice that can be at or above the freezing point. This dual representation is crucial for accurately modelling ice flow and melting processes. The model incorporates the mechanics of ice flow, including the influences of gravity, basal sliding, and internal deformation, allowing to simulate how ice responds to changes in climate, topography, and other environmental factors and making it able to be used at different scales and for different purposes (Zhang and others, Reference Zhang2013; Wang and others, Reference Wang2018).
We set the ice parameters according to Wang and others (Reference Wang, Zhang, Xiao, Ren and Wang2020), fixing the shape of the glacier and the position of the cold/temperate-ice transition zone (CTZ) according to what was interpreted with the GPR. Considering that the two glaciers are located in the same area and at approximately the same altitude, with very small variations, we set the environmental parameters of geothermal flux (56 mWm−2; according to Colgan and others, Reference Colgan, Wansing, Mankoff, Lösing, Hopper, Louden, Ebbing, Christiansen, Ingeman-Nielsen, Liljedahl, MacGregor, Á, Bernstein, Karlsson, Fuchs, Hartikainen, Liakka, Fausto, Dahl-Jensen, Bjørk, Naslund, Mørk, Martos, Balling, Funck, Kjeldsen, Petersen, Gregersen, Dam, Nielsen, Khan and Løkkegaard2022) and MAAT (according to local temperature measurements provided by a new AWS station and by Cappelen and Jensen, Reference Cappelen and Jensen2021) as the same for the two glaciers. To assess the effects of the model run time (parameter t in PoLIM), we considered different possible scenarios, but within the time range of 1–60 years, the results are indeed very similar for both WG (Figure 1S) and EG (Figure 2S). A synthesis of the data sources for PoLIM parameters is provided in Supplementary Table 2S.
4. Results
Interpretation of the GPR data sets on the two Aqqutikitsoq glaciers allowed for the estimation of their total thickness (i.e. from the topographic surface to the bedrock) at 36.35% and 47.17% of the total length of the profiles, i.e. 6.20 km and 11.93 km for the WG and EG, respectively (Figure 3S). The analysis performed by Marcer and others (Reference Marcer2017) on the western dataset (2014) yielded a higher percentage (45.04%), but here we did not interpret both where the bedrock exceeded the maximum penetration of the EM waves, and where the scattering obscured the basal reflector (i.e. where the signal-to-noise ratio was too small). In addition to the total thickness of the ice, we interpreted the EM nearly transparent zone as cold ice (CI), while the highly scattered portion as warm ice (WI), (Figure 2). To better constrain such an interpretation, we exploited some signal attributes. In particular, the Sweetness attribute (Figure 2b,f) shows more clearly than the signal amplitude the abrupt transition between CI and WI (Figure 2a,e); the Dominant Frequency (Figure 2c,g) shows overall higher values for CI when compared to WI; the Trace Envelope (a.k.a. Instantaneous Amplitude) is helpful to detect the bedrock also below high scattering zones (Figure 2d,h).
Figure 2. GPR attribute analysis. Two exemplary longitudinal profiles of the western (WG) and eastern (EG) glaciers in amplitude (a,e) are compared with three different signal attributes: sweetness (b,f); dominant frequency (c,g); trace envelope (d,h). CI cold ice; WI warm ice; l layered ice. The black line marks the glacier bottom. Vertical exaggeration 3x.
4.1. Ice thickness/volume
Figure 3a and d shows the total ice (TI) thickness along the GPR profiles obtained by interpolation in the zones where the bedrock was not interpreted and the thickness of cold (shallower) and warm (deeper) ice (CI and WI from hereafter) for both the WG and EG (Figures 3b,e and c,f, respectively). Due to the overall high quality of the data, it was possible to interpret the thickness of CI in all profiles, while we only detected WI in the central part of both glaciers, as is expected for polythermal glaciers (Glasser, Reference Glasser, Singh, Singh and Haritashya2011).
Figure 3. Ice thickness of western (WG, 2014) and eastern (EG, 2024) Aqqutikitsoq glaciers along GPR profiles. (a,d) TI thickness; (b,e) CI thickness; (c,f) WI thickness. Red lines conventionally limit the extension of the glaciers (see Figure 1).
All the data obtained from GPR are summarised in Table 2. The area of the EG (i.e. 6.43 km2) is more than double that of the WG (i.e. 2.94 km2), as it is the estimated volume of ice (3.78 vs 1.46 108 m3). The calculated volume of the EG is indeed very similar to the estimate (1.49 108 m3) provided by (Marcer and others, Reference Marcer2017) using the same 2014 GPR data set here re-processed and re-interpreted. However, the maximum ice thickness of the WG is slightly higher than that of the EG and the maximum thickness of the CI and WI for both glaciers varies considerably being 124 m and 163 m, and 140 and 153 m, respectively. This demonstrates the complex thermal regime and heat transfer mechanisms occurring into the glaciers.
Table 2. Synthesis of estimates from GPR for WG and EG datasets
Checking the thickness of the different types of ice and the TI at the intersections of the GPR profile, we see a very good agreement, with only local differences up to 9.5 m for the WG and 6.2 m for the EG.
Figure 4 shows the results obtained by interpolating the previously described thicknesses of TI and the boundary of the glaciers, where the ice thickness is, by definition, zero (Figure 4a,d). The same was done for CI (Figure 4b,e). WI was interpolated only where it was detected, i.e. only in the central part of both glaciers (Figure 4c,f). We used the Kriging algorithm with a circular search space of variable radius, automatically estimated and inversely proportional to the data density. The histogram distribution of TI thickness classes is provided in Figure 4S.
Figure 4. Interpolated ice Thickness of western (2014, WG) and eastern (2024, EG) Aqqutikitsoq glaciers from GPR profiles. (a,d) total ice thickness; (b,e) cold ice thickness; (c,f) warm ice thickness (interpolated only where present). Thickness is always set to zero at the glacier margins. Red lines conventionally limit the extension of the glaciers.
The interpolated ice thickness has the maximum values along the central part of the WG and aligned in an east-west direction or localised in the northern part for the EG. However, the spatial distribution of the GPR profiles is not homogeneous. In the WG, they are concentrated in the northern part of the glacier, while in the EG they are more homogeneously distributed, except in the peripheral portions. Instead of performing a statistical analysis of bedrock uncertainties (e.g. Marcer and others, Reference Marcer2017), which is somewhat subjective as it involves several parameters, we considered the interpolation error simply as a function of the distance to the measured profile points, following Gillespie and others (Reference Gillespie, Yde, Andresen and Citterio2023) approach (Figure 5). For the WG and EG, 16% and 31% of the areas are more than 100 m away from the measured data, respectively, but while in the former case almost all these zones are in the south, in the latter they are evenly distributed over the entire glacier surface.
Figure 5. GPR uncertainty analysis. (a,c) Distance from GPR profiles and glacier boundaries for WG and EG. (b,d) Histograms of distances of cells from data.
We compared the ice thickness results obtained from GPR (and DEM) data (Figures 4a,d and 6c,g) with the values reported in Millan and others (Reference Millan, Mouginot, Rabatel and Morlighem2022) for all the glaciers in the RGI repository (RGI, 2023), -hereafter simply referred to as RGI- (Figure 6a,e) and with the modelling obtained with the GlabTop2 module (Figure 6b,f). For both the models, we used the limits provided in the database, which are slightly different from those used to interpolate the GPR data. The mean ice thicknesses estimated for the WG and the EG are 33.9 m and 54.2 m for RGI and 65.4 m and 71.5 m with GlabTop2, while they reach 50 m and 59 m using the GPR data. Such differences are not only related to apparent differences in ice thickness variations within the glaciers, but also, at least in part, to different behaviour at the glacier margin, which by definition always has zero ice thickness.
Figure 6. Comparison between ice thickness estimates obtained using the RGI data (a,e); the GlabTop2 model (b,f) and the interpolation of the GPR datasets (c,g). The latter two maps are made with the same data as in Figure 4a,d. In d and h the thickness values along two longitudinal profiles (dashed lines in a, b, c, e, f, g) are plotted for the RGI model (in green), the GlabTop2 model (in red), and GPR (in black). see text for details.
In any case, the RGI underestimates the maximum ice thickness for both the WG and EG, whose maximum estimated thicknesses are 79.6 m and 66.1 m, respectively. Moreover, even if the maximum thickness zone for the WG is shifted to the north, as obtained from GPR data, the lateral variations for both glaciers are not properly estimated, as the thickness is severely underestimated and too smoothed. GlabTop2 results give maximum ice thicknesses for the WG and EG of 192.9 m and 187.9 m, respectively. While the latter is very similar to the GPR estimate (188 m), the former is quite different because the maximum found by GPR equals 160 m (see also Table 2 and Figure 5S). In addition, the GlabTop2 results for the WG are in fairly good agreement with the GPR results, since in both cases an increased ice thickness is obtained in the northeastern part of the glacier, with rapid spatial variations highlighted by both methods. This is not the case for the EG: while GPR data show a general increase in ice thickness along an east-west direction in the central part of the glacier (with a local spot to the north), GlabTop2 results show unrealistic lateral variations, sometimes with maxima at the glacier margin. This is even more apparent when ice thickness is analysed along a longitudinal profile for the two glaciers (Figure 6d,h). However, in the original version of GlabTop2, this is prevented by not considering equation 1 close to the margin of the glacier (Frey, Reference Frey2014). RGI data underestimate ice thickness and provide overly smoothed estimates, while GlabTop2 results are globally similar to GPR for the WG, but not for the EG. At the glacier margins, both modelling estimates show some artefacts with unrealistic local ice thickness increases.
To further evaluate the ice thickness estimates from GPR surveys, we make a comparison between the total ice volumes from GPR data and those resulting from some global-scale power-law relationships with parameters set as suggested by different authors (Table 3). Interestingly, the global-scale estimates give values that are quite similar to the GPR ones. The worst performance for the WG overestimates and underestimates the value provided by GPR by 19.2% and 19.5%, respectively; the worst performance for the EG overestimates and underestimates the value provided by GPR by 26.5% and 7.4%, respectively. For the WG, the parameters proposed by Radić and Hock (Reference Radić and Hock2010) give the closest results to GPR (0.158 vs 0.146 km3); for the EG, the closest results are obtained with the parameters of Chen and Ohmura (Reference Chen and Ohmura1990), (0.356 vs 0.378 km3).
Table 3. Comparison between volume estimates obtained for the western and eastern glaciers (WG and EG) by GPR datasets and by global-scale estimates with parameters proposed by different authors
* Magnitude after 100 years of sustained recession.
4.2. Thermal modelling
Thermal modelling was conducted along two longitudinal GPR profiles to allow a direct comparison between the CI and WI zones detected by the geophysical data and the temperature behaviour obtained by modelling. As previously noted in both the WG and EG, the GPR data showed large variations in the occurrence and thickness of CI and WI (Figures 2,7a,8a). As expected, WI is absent near the glacier snout, while it is generally thicker at higher elevations, with zones where it is either absent (WG, Figure 7a) or only a few metres thick (EG, Figure 8a). In both glaciers there is a zone not far from the glaciers’ snout where the bedrock is deeper and the WI thickness increases (Figures 7a,8a; see also Figure 4c,f).
Figure 7. GPR (a) and thermal modelling results (b, c, d) along a longitudinal profile of the Western Glacier (2014) for different maats values. CI cold ice; WI warm ice; WP water percolation; CTZ cold/temperate-ice transition zone. The black line marks the glacier bottom. Vertical exaggeration 3×.
Figure 8. GPR (a) and thermal modelling results (b, c, d) along a longitudinal profile of the eastern glacier (2024) for different maats values. CI cold ice; WI warm ice; CTZ cold/temperate-ice transition zone. The black line marks the glacier bottom. Vertical exaggeration 3×.
Since the calculated 1991–2020 MAAT is equal to −7.8°C ± 1.4°C, we modelled the thermal behaviour for −7.8°C, −6.4°C and −9.2°C (i.e. ± σ), keeping all other parameters constant (see Methods section).
The PoLIM results for the WG (Figure 7b–d) predict the development of WI only for the upper half of the glaciers, with an overall thickness comparable to that detected by the GPR data, but never able to match the high spatial thickness variability. The model closest to the GPR results seems to be the one for the warmer MAAT (Figure 7b), but also in this case the model is not able to predict the presence of WI at lower elevations, even if a local warming is modelled at a MAAT of −6.4°C, where the thicker WI is observed on GPR towards the WG snout.
The PoLIM results for the EG (Figure 7b–d) show an overall better agreement with the GPR. In the upper part of the glacier, all models for the different MAATs tested predict the development of WI, with the model calculated for −6.4°C (Figure 8b) being the most similar, although with slightly lower predicted thicknesses. Towards the snout, only the warmer model correctly predicts WI formation in the zone where it is imaged by the GPR data, but again its modelled thickness is lower.
5. Discussion
5.1. Thickness/volume modelling
We see a general good agreement between the result based on Macheret and others (Reference Macheret, Cherkasov and Bobrova1988) and Bahr and others (Reference Bahr1997) formulas and our GPR estimates, with some minor differences due to the different coefficients proposed by the different authors (Table 3). In fact, if we consider the work of Adhikari and Marshall (Reference Adhikari and Marshall2012), they argue about the different parameters to consider for the variation of (c, γ) and try to give the best coefficients according to some environmental parameters that should be known (e.g. accumulation zone, bedrock slope, valley shape, climatic parameters). The ranges of c and γ for individual glaciers are likely to be even larger than previously thought, and the consequences of glacier-specific conditions and significant climatic disequilibrium need to be better understood in order to better estimate glacier volumes from volume-area relationships at local and regional scales (Yde and others, Reference Yde2014). In any case, volume model results appear to be good for global estimates, but not for defining details of bedrock morphology and hence local ice thickness variations, which are essential for thermal modelling and understanding heat flows due to ice-bedrock interactions. There are several examples where there is generally good agreement between ice thickness estimates from GPR surveys and modelling, but larger discrepancies occur where the bedrock morphology is not smoothed but has relevant local variations (Sattar and others, Reference Sattar, Goswami, Kulkarni and Das2019; Vergnano and others, Reference Vergnano, Franco and Godio2024). In the latter article the authors used the ice-thickness modelling algorithms to retrieve the bedrock topography where it could not be detected by the GPR. This could be an interesting approach, at least in some cases. After validation of the model used where the GPR data provide the ice thickness with a high degree of accuracy, modelling results can be used where the bedrock cannot be detected by GPR or it can only be interpreted tentatively. Vergnano and others (Reference Vergnano, Franco and Godio2024) pointed out that in the presence of scattering (due to various possible causes), the use of ground-based low-frequency antennas instead of helicopter-borne relatively high-frequency antennas cannot significantly improve the detection of bedrock, so that an integrated geophysical and modelling approach may be advisable. In our study area, the modelling results, while providing quite good estimates on a global scale, are never accurate enough locally, which prevents their extensive use where the bedrock is not detected by GPR. We already pointed out that the ice thickness estimates obtained with GlabTop are overall quite similar to the GPR results, even if the modelling produces local unrealistic outliers and lateral changes, especially for the EG.
As suggested by Švinka and others (Reference Švinka, Karušs and Lamsters2025), we calculated the shape factor f (see the Methods section) along several transects of the two glaciers, then averaging the results (Table 1S). We obtained mean f values equal to 0.92 and 0.94 for the WG and EG, respectively. By applying such values instead of 0.8, which is considered the default for valley glaciers, we obtained a maximum ice thickness closer to GPR (Figure 7S) for the WG. In fact, the maximum thickness of the WG reaches 163.9 m and 187.9 m for f equal to 0.92 and 0.8, respectively, while is 159.8 m for GPR data. On the contrary, the maximum thickness of the EG reaches 165.3 m and 192.9 m for f equal to 0.94 and 0.8, respectively, while it is 188.3 m for the GPR data, so the use of the specifically calculated f factor seems to underestimate the maximum ice thickness (Table 3S).
More interestingly, the mean ice thickness values are always closer to the GPR estimates when calculated with the specifically calculated f factors for both the WG and EG. In fact, the mean ice thickness of the WG reaches 56.9 m and 65.4 m for f equal to 0.92 and 0.8, respectively, while it is 49.8 m for the GPR data. Similarly, the mean thickness of the EG reaches 60.9 m and 71.5 m for f equal to 0.94 and 0.8, respectively, while it is 58.6 m for the GPR data (Table 3S).
In addition, the adjusted f factors allow to obtain better results close to the boundaries of both glaciers, also smoothing some unrealistically high thicknesses estimated, for example, at the northwestern boundary of the WG (see Figure 7S for details).
Furthermore, in all the estimates it should be taken into account that glaciers in general do not consist only of pure ice, but they can contain a small fraction of water, air, as well as other impurities (Colucci and others, Reference Colucci, Forte, Boccali, Dossi, Lanza, Pipan and Guglielmin2015; Zhao and others, Reference Zhao, Forte, Colucci and Pipan2016), which have to be considered when calculating the volume and, more importantly, the water stored in it (i.e. the water equivalent). However, for the two study glaciers, no units other than cold and warm ice are relevant.
From the signal analysis, a remarkable effect in terms of signal penetration and signal-to-noise ratio seems to be more related to the environmental and ice conditions during the survey. In fact, the 2014 survey was conducted in August, while the 2024 in February, with an average temperature about 25°C lower and certainly less liquid water at the glacier surface, which favoured better GPR performance.
Regarding the spatial resolution of the GPR datasets, both have limited data south of the glaciers due to logistical and safety constraints. This lack of data certainly affects the interpolation with the glacier margin, most likely resulting in an underestimation of ice thickness near the southern margin. The use of helicopter-borne GPR equipment (e.g. Rutishauser and others, Reference Rutishauser, Maurer and Bauder2016; Santin and others, Reference Santin2019) or unmanned aerial vehicles (UAVs) (e.g. Ruols and others, Reference Ruols, Baron and Irving2023) could partially overcome this problem, allowing data to be collected even in some additional areas. However, even airborne instruments have some inherent limitations, and special care must be taken when collecting and processing data (Forte and others, Reference Forte, Basso Bondini, Bortoletto, Dossi and Colucci2019).
Other geophysical methods can contribute to the detection of the glaciers bed, both when the ice thickness is too high for GPR and when the warm ice limits the penetration of EM waves. Active (Navarro and others, Reference Navarro, Macheret and Benjumea2005; Johansen and others, Reference Johansen2011) and passive (Picotti and others, Reference Picotti, Francese, Giorgi, Pettenati and Carcione2017) seismic methods have been shown to be applicable on glaciers, despite some inherent logistical and resolution constraints. Recently, Distributed Acoustic Sensing (DAS) has been successfully tested for glaciological investigations exploiting surface (Manos and others, Reference Manos2024) and borehole (Booth and others, Reference Booth2023) measurements, showing its potential for possible different future investigations.
5.2. Thermal modelling
The results of the thermal modelling are only partially consistent with the WI and CI distributions as obtained by the GPR data. It is well known that the thermal structure of ice bodies reflects both environmental setting and internal processes (Blatter and Hutter, Reference Blatter and Hutter1991). The thermal models typically do not consider geophysical geometric constraints, with some exceptions. Wilson and others (Reference Wilson, Flowers and Mingo2013) used GPR and borehole temperature data for two valley glaciers in the Yukon, Canada coupled with thermal modelling results. Delcourt and others (Reference Delcourt, Van Liefferinge, Nolan and Pattyn2013) mapped the CTS on GPR data, validating the results with borehole temperatures and 2D thermal models simulating different glacier conditions for the McCall Glacier, Alaska, USA with the goal of defining the past and future evolution of the Equilibrium-Line Altitude (ELA). However, in both cases the thermal modelling is performed on profiles with limited spatial variability in ice thickness and only along longitudinal profiles of valley glaciers with very regular bedrock geometry. In particular, Wilson and others (Reference Wilson, Flowers and Mingo2013) concluded that the stronger thermal effects are related to the meltwater trapped in the study glaciers, while the strain heating generally plays a minimal role. This is not the case for the two studied Aqqutikitsoq glaciers, where the spatial variability of the WI and CI distribution is most likely related to local higher and lower strain conditions, in turn due to the rough topography of the bedrock and to the very complex shape of the glaciers (see map in Figure 1). Furthermore, the very high spatial thermal variability of the two study glaciers is indeed not imaged by GPR datasets collected on other similar glaciers both in E-Greenland (Citterio and others, Reference Citterio, Mottram, Hillerup Larsen and Ahlstrøm2009; Abermann and others, Reference Abermann, Van As, Petersen and Nauta2014) and in W-Greenland (Gillespie and others, Reference Gillespie, Yde, Andresen and Citterio2023), as well as on glaciers at similar latitudes in Iceland (e.g. Lamsters and others, Reference Lamsters, Karušs, Rečs and Bērziņš2016) and in the Svalbard Islands (see e.g. the comprehensive dataset reported in Sevestre and others, Reference Sevestre, Benn, Hulton and Bælum2015). This behaviour may be explained by the highly irregular bedrock of the Aqqutikitsoq glaciers, or by a possible ongoing thermal switch driven by different variables and conditions, as discussed in detail in Sevestre and others, Reference Sevestre, Benn, Hulton and Bælum2015.
Interestingly, when we set the upper (i.e. warmer) limit of MAAT variability in PoLIM, we obtained modelling results that were more compatible with the WI and CI distributions detected by GPR (Figures 7, 8). Even considering the high variability of the MAAT in the 30-year period (1991–2020), with the minimum value of −10.5°C (1992) and the maximum of −3.8°C (2010), a warming trend is apparent (Figure 8S). We could therefore hypothesise that what is recorded by GPR for the WG (August 2014) and EG (February 2024) images warmer conditions than the 30-year mean MAAT. Due to the irregular shape and bedrock of the study glaciers and the lack of borehole temperature time series, it is not possible to validate this hypothesis and make consistent quantitative predictions for the ELA and thermal evolution of the glaciers. Furthermore, we should certainly quantify the feedback of the avalanches irregularly feeding the WG and EG, and we cannot straightforwardly predict if these two glaciers will have a transition from polythermal to completely cold (in a counter-intuitive way) before disappearing (Delcourt and others, Reference Delcourt, Van Liefferinge, Nolan and Pattyn2013; Wilson and others, Reference Wilson, Flowers and Mingo2013) under ongoing warming climate and the overall long-term negative mass balance of observed on western Greenland glaciers (Securo and others, Reference Securo2024).
6. Conclusions
While global-scale ice thickness models provide relatively accurate estimates of regional glaciers volume, they have significant limitations in local accuracy, particularly where the bedrock exhibits high spatial variability. This is evidenced where the models are producing local outliers and unrealistic changes, highlighting the challenges of consistently applying these models in areas of complex glacier topography.
We found that integrated GPR signal attributes can be very useful to better image both CI and WI and to better detect the bedrock, even below the high scattering. In fact, attributes can better detect different EM facies that represent the cumulative signature of a material on GPR data. Very recent and innovative attribute applications to sedimentary environments (Koyan and Tronicke, Reference Koyan and Tronicke2024) and using deep learning algorithms (Roncoroni and others, Reference Roncoroni, Forte and Pipan2024) could be conveniently exported to glaciological datasets.
Our study found that 2D thermal modelling results were only partially consistent with WI and CI distributions derived from GPR data, indicating the unique and complex thermal structures of polythermal glaciers with complex shape and geometry. This reflects the not obvious balance between environmental influences and internal processes. Thermal models often overlook critical geophysical constraints, resulting in an incomplete representation of the diverse glacier conditions observed in our study. Furthermore, 2D modelling obviously cannot handle the 3D variability, which may lead to oversimplified conclusions.
Our results also highlight the critical role of geophysical field measurements in constraining bedrock topography to improve model accuracy and reduce uncertainty. In general, the accuracy of current ice thickness models is highly dependent on the availability of in situ data (Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022), which remains inadequate for most glaciers (Welty and others, Reference Welty2020).
In addition, the observed 0.6°C/decade warming trend between 1991 and 2023 in the study area suggests that the thermal conditions recorded by GPR in certain glaciers may reflect a deviation from the historical mean. However, given the irregular shape of the study glaciers and the lack of borehole temperature records, it remains difficult to validate this hypothesis or to make accurate predictions regarding the Equilibrium Line Altitude (ELA) and long-term thermal evolution of these glaciers.
Our findings highlight the need for improved regional-scale inversion methods and high-resolution field data to better capture the spatial variability in glacier flow, geometry and mass balance, especially in small and/or complex glaciers. Filling these data gaps is in turn essential to accurately reconstruct the temporal evolution of global ice reservoirs and to understand the impact of climate change on glacier stability and potential thermal transitions and geometrical changes.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10067.
Acknowledgements
This study is part of the Local Glaciers Sisimiut Project (LOGS), funded by a Research Promotion Grant from the Greenland Research Council (Nunatsinni Ilisimatusarnermik Siunnersuisoqatigiit). The authors would like to thank the Greenland Winter Warning Association (GWWA) for their support during fieldwork, particularly Arne Hardenberg, Krister Støvlbæk and Miki Inuk Reimer Geisler. The authors would like to thank the Arctic DTU and the Geological Survey of Denmark and Greenland (GEUS) for their help in the logistical support. Esplora Srl is thanked for facilitating the 2024 field data acquisition. Schlumberger is acknowledged for the Petrel academic grant provided to the University of Trieste.
Open access
