Ice thickness distribution of all Swiss glaciers based on extended ground-penetrating radar data and glaciological modeling

Abstract

Accurate knowledge of the ice thickness distribution and glacier bed topography is essential for predicting dynamic glacier changes and the future developments of downstream hydrology, which are impacting the energy sector, tourism industry and natural hazard management. Using AIR-ETH, a new helicopter-borne ground-penetrating radar (GPR) platform, we measured the ice thickness of all large and most medium-sized glaciers in the Swiss Alps during the years 2016–20. Most of these had either never or only partially been surveyed before. With this new dataset, 251 glaciers – making up 81% of the glacierized area – are now covered by GPR surveys. For obtaining a comprehensive estimate of the overall glacier ice volume, ice thickness distribution and glacier bed topography, we combined this large amount of data with two independent modeling algorithms. This resulted in new maps of the glacier bed topography with unprecedented accuracy. The total glacier volume in the Swiss Alps was determined to be 58.7 ± 2.5 km³ in the year 2016. By projecting these results based on mass-balance data, we estimated a total ice volume of 52.9 ± 2.7 km³ for the year 2020. Data and modeling results are accessible in the form of the SwissGlacierThickness-R2020 data package.

1. Introduction

Since the end of Little Ice Age around 1850, more than half of the glacier volume in the Swiss Alps has been lost, with increasing rates in the recent decades (Zemp and others, 2015). Projections of future glacier change imply that under high-emission scenarios, the Alps will be largely deglacierized by 2100 (Zekollari and others, 2019). This will lead to remarkable changes in the annual regime of river runoff (Braun and others, 2000; Milner and others, 2009; Huss, 2011). Since glaciers act as a water reservoir at diverse timescales (Beniston and others, 2018), they are important, for example for electricity production from hydropower (Beniston, 2012; Finger and others, 2012; Gaudard and others, 2016), which is relevant for Switzerland as hydropower contributes with around 55% to the total annual electricity production (BFE, 2019). Glaciers are also of high societal relevance in many other fields, such as agriculture (Fink and others, 2004; Beniston, 2012), natural hazards (Huggel and others, 2004; Haeberli, 2005), tourism (Fischer and others, 2011; Beniston, 2012), supply of drinking water (Wyer, 2008) or with regard to sea level rise (Zemp and others, 2015). On the contrary, new chances might arise, for example for the electricity production from hydropower, due to new lakes forming in deglacierizing regions (Haeberli and others, 2016; Farinotti and others, 2019), where potentially new storage power plants can be installed (Ehrbar and others, 2018).

Accurate knowledge about today's ice volume and the subglacial topography is key to calculate future changes in river runoff regimes (e.g. Gabbi and others, 2012; Gaudard, 2015) or predicting the occurrence of alpine lakes forming after glaciers have retreated (Linsbauer and others, 2012; Haeberli and others, 2016). At spatial scales beyond those of individual glaciers, this information is typically obtained from ice thickness models or scaling approaches (Bahr and others, 2015; Farinotti and others, 2017). Various ice thickness models have recently been developed which allow estimating the ice thickness distribution from the glacier surface topography, glacier outlines and sometimes auxiliary data such as mass-balance estimates, ice flow velocities or branch lines. Currently applied models were, for example, presented by Farinotti and others (2009b, extended later to the global scale by Huss and Farinotti, 2012), Morlighem and others (2011), Paul and Linsbauer (2012), Clarke and others (2013), Brinkerhoff and others (2016) and Rabatel and others (2018). Based on such models, ice thickness and volume estimates specifically for Switzerland were presented by Farinotti and others (2009a), Linsbauer and others (2012) and Fischer and others (2015), whereas additional independent regional estimates at the scale of the European Alps are available from Haeberli and Hölzle (1995) and Farinotti and others (2019).

Ice thickness estimates obtained from modeling exhibit substantial uncertainties unless they are calibrated to ground-truth data (Farinotti and others, 2021). For the Antarctic and Greenlandic ice sheets, measured ice thickness data have been acquired systematically with a comprehensive coverage (e.g. Dowdeswell and Evans, 2004; Allen and others, 2010), enabling to interpolate the ice thickness distribution in these regions more accurately (e.g. Huss and Farinotti, 2014; Morlighem and others, 2020). Andreassen and others (2015) used measured ice thickness for 32% of the glacier area of mainland Norway, and Helfricht and others (2019) relied on direct thickness observations for 46% of Austria's glacierized areas to obtain complete ice thickness maps in combination with calibrated ice thickness models. A similar study has been accomplished by Pelto and others (2020) for the Canadian part of the Columbia River Basin. For other regions, available ice thickness observations are usually too sparse and models of entire mountain ranges are only calibrated for a few glaciers and the calibration factors are then projected to all other glaciers (Farinotti and others, 2009a; Linsbauer and others, 2012; Werder and others, 2020). To improve ice volume estimates and to obtain the subglacial topography with higher accuracy, it is therefore substantial to provide more direct measurements of the glacier thickness as for instance concluded by Farinotti and others (2017).

Apart from drilling boreholes to the bedrock, which is only feasible in limited numbers, radar sounding has been established as the most common technique to measure the ice thickness of glaciers (Plewes and Hubbard, 2001; Woodward and Burke, 2007). This is true for both, environments with non-temperate ice, such as in Antarctica or Greenland (e.g. Bingham and Siegert, 2007), as well as in regions with temperate glaciers such as in Alaska, Patagonia, the Himalaya or in the European Alps (e.g. Blindow and others, 2012). For surveying temperate glaciers with radar, a major challenge arises from the substantially higher dampening of electromagnetic waves in comparison with cold ice (Watts and England, 1976). For this reason, impulse radar systems have been employed, hereafter referred to as ground-penetrating radars (GPR). They are able to transmit ~1 cycle at high voltages and with a relatively low central frequency defined by the antenna length (Plewes and Hubbard, 2001). Typically, their frequencies range from a few megahertz, depending on the maximum manageable antenna length, up to several tens of megahertz, depending on the required maximum penetration depth.

Airborne GPR systems with sufficiently low frequency to penetrate temperate glaciers have been used worldwide since decades (Watts and Wright, 1981; Kennett and others, 1993; Christensen and others, 2000; Conway and others, 2009; Blindow and others, 2011; Rignot and others, 2013; Gacitúa and others, 2015; Urbini and others, 2017; Langhammer and others, 2019b; Pritchard and others, 2020), carried either by fixed-wing airplanes on wide and flat glaciers and ice caps or by helicopters on glaciers in narrow valleys or high-altitude environments with steep glaciers. Also in the Swiss Alps, several helicopter-borne GPR systems have been employed during the last two decades as summarized by Rutishauser and others (2016). Together with a number of ground-based surveys (e.g. Bauder and others, 2003; Huss and Fischer, 2016; Feiger and others, 2018), they contribute largely to the ice thickness datasets for this region, and are available through the public glacier ice thickness database GlaThiDa (Gärtner-Roer and others, 2014; Welty and others, 2020). To complement this pre-existing database, we have established the new helicopter-borne GPR platform AIR-ETH (Langhammer and others, 2019b). It is the first GPR platform with two orthogonal antenna pairs, which significantly improves the bedrock recovery of valley glaciers (Langhammer and others, 2017).

Our study has three objectives: first, we present a large amount of new ice thickness data from Swiss glaciers, acquired using the AIR-ETH system. These data substantially extend the already existing database of ice thickness point measurements. Second, we generate continuous maps of the ice thickness distribution through interpolation by glaciological modeling based on the extended and highly comprehensive set of ice thickness data. For that purpose, we apply two different algorithms, namely the recently developed Glacier Thickness Estimate (GlaTE) algorithm (Langhammer and others, 2019a) and the Ice Thickness and Volume Estimation method Observation-based (ITVEO) that is based on the concepts of Huss and Farinotti (2012). Additionally, we use the ice thickness distributions and the corresponding GPR data to calibrate the models for estimating the ice thickness distribution of the remaining, mostly small glaciers for which still no GPR data exist. This is the first time that at the regional scale, ice thickness distribution and glacier bed topography has been determined based on such a high density of direct ice thickness observations. The third and final aim of our study is to provide a new estimate of the total ice volume in the Swiss Alps.

2. Study region and datasets

A glacier inventory is a prerequisite for a regional assessment of ice thickness distribution and glacier volume. Here, we refer to the Swiss Glacier Inventory 2016 (SGI2016, GLAMOS, 2020b; Linsbauer and others, subm.) for the definition of glacierized areas. This inventory comprises 1400 glaciers with a total area of 961 km² and refers to the years between 2013 and 2018, with a center date in 2016. A number of 160 of these glaciers are larger than 1 km² and cover a total area of 767 km², whereas the overwhelming majority of the glaciers are small and, thus, contribute little to the total ice volume (e.g. Fischer and others, 2014). Most glaciers are located in the Bernese and Valais Alps, where the Alps reach the highest elevation within our study region. Substantial glacierized areas also exist in the central and eastern Swiss Alps. The outlines provided by the SGI2016 define the lateral extent of our glacier maps. Additionally, the SGI2016 provides raster data indicating which parts of the glaciers are covered by supraglacial debris (GLAMOS, 2020b).

Our study region is covered by swissALTI3D, a recent, high-quality digital elevation model (DEM) provided by the Federal Office of Topography (Swisstopo, 2010). It is accompanied by well-documented meta data. The DEM used for modeling is the 2019-release of the swissALTI3D with 2 m-resolution (Swisstopo, 2019), which we down-sampled to a resolution of 10 m. The acquisition year over the extent of all glaciers can be retrieved from the meta data, which is provided within cells of 1 km × 1 km. It provides the mode (most occurring) year, as well as the youngest and oldest recording years of the base data used for producing the DEM (Weidmann and others, 2018). On average, the 2019-release of swissALTI3D represents the surface topography of the year 2015, with oldest being from 2007 and youngest from 2017.

To date, a large amount of ice thickness data for Swiss glaciers has been acquired with GPR by various researchers. Most of them are available through GlaThiDa (GlaThiDa Consortium, 2019) operated by the World Glacier Monitoring Service (Gärtner-Roer and others, 2014; Welty and others, 2020). In addition to these data, we also use the ice thickness data measured with seismic methods by Thyssen and Ahmad (1969), covering an area of Konkordiaplatz where the three main branches of Aletschgletscher merge. At this location, the largest ice thickness in the Swiss Alps occur, but to date the glacier bed underneath the deepest ice has not been detected with GPR. The overdeepening detected by Thyssen and Ahmad (1969) was later confirmed by drilling boreholes as described in the study of Hock and others (1999). In total, the data published through GlaThiDa comprise around 950 km of GPR profiles, mostly acquired on the larger glaciers in the Bernese and Valais Alps. Additionally, the Glaciology Group at the Laboratory of Hydraulics, Hydrology and Glaciology (VAW, ETH Zurich) archived another ~550 km of GPR profiles, which were measured by various researchers but not published so far. References of the corresponding studies are listed in Table 2. The ice thickness data were complemented with another 950 km of profiles, recorded in the framework of this study as introduced below. They were recorded not only in the Western Swiss Alps, where most of the previous surveys were conducted, but also on the glaciers in the central and eastern Swiss Alps, for which almost no data existed before. An overview of which glaciers have been surveyed to date is shown in Figure 1a. In Figure 1b, the location of GPR sections is shown for the Bernese Alps. For analogous figures for other regions, refer to the Supplementary material.

Fig. 1. Overview of GPR coverage: (a) SGI2016 outlines of all glaciers in the Swiss Alps. Glaciers for which no GPR data have been acquired to date are in blue. Glaciers for which GPR data are available are shown in red, with the color-code indicating the density of the data coverage expressed by the mean of the distances of all points on the glacier to the closest GPR-measurement point. (b) Zoom into the Bernese Alps with the location of the GPR-measurement points shown in light gray to black depending on the data source. See the Supplementary material for other parts of the Swiss Alps. Background: hillshade from the swissALTI3D after Swisstopo (2019).

Since ice thicknesses were measured in various years and seasons, ice thicknesses used as input for the ice thickness distribution modeling were standardized by subtracting the measured glacier bed from the DEM. This requires the elevation datum of measured glacier surface- and glacier bed-elevation to be consistent with the datum of the DEM. Therefore, we compared the surface elevation accompanying the GPR datasets with a DEM dating back as close as possible to the time of GPR-data acquisition. In case the two surface elevations differed strongly, which occurred when surface elevations were not measured simultaneously during GPR acquisition, directly the measured ice thickness was used as input for the ice thickness interpolation. This results in increased uncertainties for the interpolated ice thickness, which was taken into account as detailed in Appendixes C and D. This was the case for about a dozen small glaciers surveyed between 2009 and 2015 (see metadata of SwissGlacierThickness-R2020, Grab and others, 2020).

3. Methods

3.1. GPR surveying

3.1.1. GPR data acquisition

GPR data were acquired in the winter and spring seasons of the years 2016–20 with the GPR platform AIR-ETH (Langhammer and others, 2019b). This GPR platform consists of the pulseEkko radar system from Sensors & Software and is carried by a helicopter as a sling-load. For the campaigns presented here, it was equipped with two orthogonal transmitter–receiver antenna pairs of 25 MHz center frequency, in order to record dual-polarization data. Both pairs were mounted in broadside configuration and were recording alternatingly around eight traces per second each. With the target flight speed of 37 km h⁻¹ (20 knots), this resulted in a trace every 0.8 m available for the consequent data processing. From three global navigation satellite system (GNSS) antennas, the pitch, roll and yaw of the platform was obtained, whereas a laser-altimeter constantly logged the distance to the glacier surface. From this information, the glacier surface elevation underneath the platform was deduced. The recording and processing of the differential GNSS and laser altimeter-data, as well as the navigation aid for the pilot, was performed by GEOSAT SA, a Swiss-based surveying company. The target height of the antennas above ground was defined as 30 m, which in field-implementation typically varied by around ±20 m. Larger heights were avoided for data-quality reasons and lower heights for safety reasons.

For typical valley glaciers, one longitudinal (i.e. along the glacier flow direction) and several cross-profiles, around every 250 m (small glaciers) and up to every 750 m (larger glaciers), were recorded. For more complex-shaped glaciers, this configuration was adapted with a similar profile density, seeking roughly orthogonal intersections. This facilitated the discrimination of off-plane reflections and reflections from steep valley walls by comparing orthogonal profiles during data interpretation. Profile locations (interpreted parts) are shown for an example region in Figure 1b. The ultimate goal of our campaigns was to cover as many glaciers as possible with the given finances. Therefore, data acquisition missions were conducted region-wise and it was aimed to cover possibly all glaciers in each region, unless they were geographically very isolated and/or very small.

3.1.2. GPR data processing

A data-processing workflow was established based on the experience of previous studies with helicopter-borne GPR (Merz and others, 2015; Rutishauser and others, 2016; Grab and others, 2018; Langhammer and others, 2019b). It mainly consists of standard GPR-processing steps that are also commonly used for processing ground-based data. A crucial adaption, however, is required because data from helicopter-borne GPR are contaminated with severe ringing noise (Fig. 2a) due to the interference of the electromagnetic waves with the helicopter (Rutishauser and others, 2016; Langhammer and others, 2019b). All processing steps were implemented in our in-house processing package GPRglaz (e.g. Grab and others, 2018), which is based on Matlab and uses parts of the CREWES library (e.g. Margrave and Lamoureux, 2019).

Fig. 2. Selected data-processing steps for profile 7 recorded on Kanderfirn (A55b/13 in Fig. 1b) in February 2019: (a) rawdata, (b) data after ringing removal and bandpass filtering and (c) data after focusing and time-to-depth conversion using Kirchhoff migration.

Data from an example GPR-profile are shown in Figure 2 after selected processing steps. It has been recorded in February 2019 and its location is marked as Profile 7 in Figure 1b. The main processing steps are (1) a time zero correction based on the arrival of the direct wave, (2) removal of the ringing noise using an optimized procedure of singular value decomposition filtering as described in detail by Grab and others (2018), (3) bandpass filtering, (4) merging of the dual-polarization data channels as described in detail by Langhammer and others (2019b), (5) trace binning for a regular sampling at a spatial sampling rate depending on the flight speed and (6) image focusing and time-to-depth conversion by migrating the data with constant radar wave velocities for ice and air of 0.169 m ns⁻¹ (Glen and Paren, 1975) and 0.299 m ns⁻¹, respectively. For the migration, a Kirchhoff time migration scheme (after Margrave and Lamoureux, 2019) was used which yielded satisfactory results after relatively short computation times, but for some profiles also a computationally more expensive reverse time migration was considered (Grab and others, 2018). The identical processing flow was applied to data from both channels of the dual polarization AIR-ETH system. Additionally, the data of the two channels were merged by summation, and the processing was applied to the resulting trace as well. Therefore, three images were available for the interpretation: one from records with the antennas orthogonal to flight direction (channel 1), one from records with the antennas parallel to flight direction (channel 2) and one for the merged channels.

3.1.3. Interpretation of glacier bed reflections

For many profiles, images from the merged channels show substantial improvements in quality (see Langhammer and others, 2019b, for more details). In some cases, it turned out to be useful to additionally consider the individual channels during data interpretation. In a number of profiles the glacier bed reflection was only visible in the data of one channel, thereby confirming the usefulness of the dual-polarization AIR-ETH system. Another advantage of the dual-polarization data is that unwanted events in the data, such as reflections from steep valley walls, were often only visible in one channel. An example is shown in Figures 3a, b, where some reflections with a linear move-out and slightly elevated frequency content, typical for reflections from valley walls, only occurs in the image recorded with channel 1. Such events may interfere with the glacier bed reflection and hinder the interpretation.

Fig. 3. Glacier bed interpretation of profile 7 from Kanderfirn (A55b/13): migrated reflection image recorded with (a) channel 1 and (b) channel 2, (c) glacier bed interpretation and classes depending on reflection quality and (d) corresponding ice thickness and uncertainty ranges used as input for GlaTE and ITVEO.

When interpreting glacier bed reflections, we distinguished three quality classes (Fig. 3c) in order to later account for the data quality when estimating the uncertainty in the measured ice thickness and to deal with ambiguous reflections. The first and highest quality class was assigned to coherent reflections with a clear onset. The second class was assigned for reflections not occurring coherently over greater distances or with a poorly defined onset. In either case, an alternative event could be selected, for example if it was not possible to clearly distinguish between a true glacier bed reflection occurring from within the profile and an image of an off-plane reflection. The hierarchy of main and alternative reflection was defined on a case-by-case basis, either by the coherence of the reflection or from a glaciological point of view. In cases of multiple reflections, for which no clear onset can be identified, but which most probably originate from the glacier bed, it was not possible to pick an individual event. Instead, a minimum–maximum time window was picked and labeled as third-class pick, whereas the actual ice thickness was defined by the mean of the selected range.

In case of an unambiguous interpretation of glacier bed reflections, the uncertainty in measured ice thickness was estimated from the error in the GPR wave velocity in temperate ice and the uncertainty in picking reflections in the radargrams. For the GPR wave velocity we used 0.169±0.005 m ns⁻¹ following Glen and Paren (1975), which is less conservative than the uncertainty of ±0.008 m ns⁻¹ used by Lapazaran and others (2016) as that study applied it to a polythermal glacier. For the picking uncertainty, Rutishauser and others (2016) estimated an uncertainty of ±5 m. From our experience, this picking uncertainty is an appropriate estimate on average. Here, we used ±3 m for reflections with a well-defined onset (class 1) and ±7 m for reflections with a poorly defined onset (classes 2 and 3). For data from the literature or from older studies, for which we do not know the data quality, an average picking uncertainty of ±5 m was assumed. In cases, where the interpretation was ambiguous and an alternative reflection was identified or reflections were picked over a time window (class 3), the uncertainty bound spans over the entire range of the identified glacier bed reflections, additionally to the velocity and the picking uncertainty. This procedure results in uncertainty ranges $\Delta h_{{\rm GPR}}^{-}$ toward smaller and $\Delta h_{{\rm GPR}}^{ + }$ toward larger ice thickness as displayed in Figure 3d.

Since the GPR data were acquired in different years and during different seasons, the measured ice thicknesses were not directly compatible with the DEM on which the ice thickness interpolation was based on. Thanks to the differential GNSS and laser-altimeter measurements with the AIR-ETH system, however, the precise elevation datum of the glacier surface and likewise of the glacier bed is known for each measured ice thickness. It was therefore permissible to subtract the glacier bed elevation from the DEM (SwissALTI3D-r2019) and use this standardized ice thickness h _GPR (see Fig. 3d) as input for the ice thickness interpolation as explained in Section 2.

3.2. Interpolation with ice thickness modeling

To obtain continuous ice thickness maps from the partly sparse coverage of GPR profiles, we applied an interpolation using two independent distributed ice thickness modeling algorithms. The basic concept is to optimize an ice thickness model using the available data for each glacier and to use the calibrated model to optimally estimate ice thickness distributions in unmeasured regions of the glacier. To this end, we rely on the algorithm Glacier Thickness Estimate (GlaTE) that has been recently developed by Langhammer and others (2019a), and on the Ice Thickness and Volume Estimation method Observation-based (ITVEO) method, a further development of the approach presented in Huss and Farinotti (2012).

GlaTE and ITVEO are well suited for the current study, because they can efficiently account for large amounts of ice thickness measurement data and are relatively inexpensive for applications to large sets of glaciers. Both models follow a concept presented first by Farinotti and others (2009b), which bases on mass conservation, apparent mass balances and Glen's ice flow law (Glen, 1955). A main difference persists in the way this concept is implemented, especially with computations being performed within ice flow sheds for GlaTE and within elevation bands for ITVEO. Furthermore, the schemes with which ice thickness observations are incorporated are also fundamentally different among the two models as explained in detail below.

The final maps of ice thickness distribution and glacier bed topography were obtained in a last step by averaging the results of the two algorithms, and is subsequently referred to as the MEAN. Unless specified differently, we hereafter refer to the MEAN result when discussing ice thickness distributions and glacier bed topography maps.

3.2.1. Glacier thickness estimate (GlaTE)

The GlaTE algorithm used in this study is identical to the one presented in Langhammer and others (2019a), except for a few modifications. GlaTE includes a modeling algorithm equivalent to the one of Clarke and others (2013). The approach builds on a physically based ice flow model that estimates the glacier bed topography from the surface topography, apparent mass-balance fields, and a stress–strain relationship for basal shear derived from Glen's law (Glen, 1955). The resulting ice thickness model is combined with the ice thicknesses measured with GPR using an inversion scheme. As detailed in Appendix A, this scheme is an extended form of the one presented by Langhammer and others (2019a), with the main difference that it also accounts for glacier bed topography gradients along the glacier boundary.

3.2.2. Ice thickness and volume estimation – observation-based (ITVEO)

The basis of the ITVEO approach for inter- and extrapolating point ice thickness observations to the entire glacier is an ice thickness model that was originally developed for global-scale applications by Huss and Farinotti (2012). In the framework of ITMIX2 (Farinotti and others, 2021) the approach has been extended to be able to incorporate thickness observations for constraining the model solution to match measurements. This methodology has been further developed in the current study as described in Appendix B.

3.2.3. Final ice thickness distribution and glacier bed topography

The two independent models GlaTE and ITVEO result in ice thickness distributions h _glate(x,  y) and h _itveo(x,  y). Since both were computed within identical glacier outlines and produced on identical grids with a 10 m resolution, no further processing was required prior to averaging the two models. With the above notation, the final maps of ice thickness distributions and glacier bed topography is simply obtained from

(1)$$h ( x,\; y) = {h_{\rm glate} ( x,\; y) + h_{\rm itveo}( x,\; y) \over 2},\; $$

to which we refer to hereafter as MEAN. It has been concluded from the Ice Thickness Model Intercomparison eXperiments – ITMIX (Farinotti and others, 2017) and ITMIX2 (Farinotti and others, 2021), that model ensembles have a higher skill in predicting ice thickness than any individual model, provided all models perform similarly well. In ITMIX2, GlaTE and ITVEO both showed a similarly good performance. Thus, averaging these two models is a first step toward a more robust areal interpolation.

The measured ice thicknesses used as input for obtaining the ice thickness distributions were standardized to be compatible with the DEM as introduced in Section 2. Therefore, DEM, ice thickness distribution and glacier bed topography form a consistent triplet. This permits the glacier bed topography to be determined by subtracting the ice thickness from the DEM. The glacier bed topography is thus consistent with the DEM and, by combining the two, terrain models of the hypothetically completely deglacierized landscape can be generated.

3.3. Ice volumes, homogenization to a specific year and temporal extrapolation

The ice volume of individual glaciers was obtained from summation of interpolated thicknesses at all gridcells, multiplied with the cell area 10 m × 10 m. The resulting ice volumes refer to the recording years of the base data on which the glacier outlines and the surface topography given by the DEM are based on, hereafter referred to as the inventory year. In order to provide a nation-wide ice volume estimate referring to a specific year, a temporal homogenization of the glacier-specific volume to a single year was required. The center year of the 2019-release of swissALTI3D is 2015 (slightly younger according to the area-weighted average), whereas it is 2016 for the SGI2016 (GLAMOS, 2020b). We thus decided to homogenize the volume to the year 2016, referred to as the reference year hereafter.

For this homogenization we rely on both long-term geodetic ice volume changes available for each individual Swiss glacier over the period 1980–2010 (Fischer and others, 2015), and in situ observations of year-to-year mass-balance variability at up to 20 glaciers in all glacierized regions of Switzerland spanning the period 1970–2020 (GLAMOS, 1881–2018, 2019). For each glacier with direct observations, the annual mass-balance anomaly with respect to that glacier's long-term mean, given by Fischer and others (2015), was evaluated and spatially extrapolated to all unmeasured Swiss glaciers based on an inverse-distance weighting scheme. For extrapolation to glaciers North of the main Alpine crest this scheme attributes a higher weight to observed mass balances on the same side of the weather divide, and vice versa for glaciers south of the main Alpine crest. Annual mass-balance anomalies ΔB _i,y of glacier i and year y were then superimposed on the glacier's long-term average mass balance $\overline {B}_i$. Converting glacier mass balance into an annual change in ice volume requires an estimate of the time-evolving area A _y for each glacier and a mean density ρ_ΔV = 850 kg m⁻³ (Huss, 2013) for the occurring volume change. Mean annual glacier-specific area-change was derived from the observed glacier area change between 1973 and 2010 (Fischer and others, 2014). The same rate was also used to extrapolate glacier area up to 2020. Ice volume change ΔV _g,y of glacier i in year y is thus obtained as

(2)$$\Delta V_{i, y} = ( \Delta B_{i, y} + \overline{B_i}) \cdot A_{i, y} \cdot {\rho_{\rm {water}}\over \rho_{\rm {\Delta V}}}.$$

Summation of annual ice volume changes relative to ice volume for the inventory year, that is the volume based on direct thickness observations, yields series for the evolution of each glacier's volume over time. Overall, ice volume in Switzerland is finally obtained for every year between 1970 and 2020 as the total volume of all glaciers. This series is useful to relate our new regional ice volume estimate to previous studies, and to analyze relative changes in ice volume.

4. Uncertainties

4.1. Point-specific ice thickness uncertainty

The point-specific ice thickness uncertainty is estimated as the standard uncertainty u ^±(x,  y) in ice thickness at an arbitrary point of the resulting ice thickness models. For its estimation, we consider uncertainties of the ice thickness interpolation $u_{{\rm int}}^{\pm }( x,\; \, y)$, uncertainties in GPR measurements $u_{{\rm gpr}}^{\pm }( x,\; \, y)$ and uncertainties in surface elevation $u_{{\rm surf}}^{\pm }( x,\; \, y)$. The superscript ‘+’ indicates uncertainty toward larger and the superscript ‘−’ toward smaller ice thickness, respectively. Quantity $u_{{\rm int}}^{\pm }( x,\; \, y)$ was estimated using a bootstrap technique (Efron, 1979), $u_{{\rm gpr}}^{\pm }( x,\; \, y)$ by repeating the ice thickness interpolation while replacing the measured ice thickness input h _GPR by $h_{{\rm GPR}}-\Delta h_{{\rm GPR}}^{-}$ and h _GPR+$\Delta h_{{\rm GPR}}^{ + }$, and for $u_{{\rm surf}}^{\pm }( x,\; \, y)$ we estimated a value which is constant across individual glaciers. More details about the computation of these uncertainty components are given in Appendix C.

Since quantities $u_{{\rm int}}^{\pm }( x,\; \, y)$, $u_{{\rm gpr}}^{\pm }( x,\; \, y)$ and $u_{{\rm surf}}^{\pm }( x,\; \, y)$ are controlled by mechanisms independent from each other, they are presumed to be uncorrelated. The total ice thickness uncertainty at a specific point is therefore obtained from

(3)$$u^{\pm}( x,\; y) = \sqrt{u_{{\rm int}}^{\pm}( x,\; y) ^{2} + u_{{\rm gpr}}^{\pm}( x,\; y) ^{2} + u_{{\rm surf}}^{\pm}( x,\; y) ^{2}},\; $$

which is applied separately for quantities u ⁺ and u ⁻.

4.2. Ice volume uncertainty

The ice volume uncertainty $\bar u_{{\rm V}}$ is expected to be substantially smaller than the sum of the point-specific ice thickness uncertainties u ^±(x,  y) introduced above, because values of u ^±(x,  y) cancel each other out to a certain extent when averaging over entire glaciers and, even more, over the entire Swiss Alps. Since uncertainties in total glacier volumes can be regarded as mean uncertainties, we use the overline-notation to distinguish the volumetric uncertainty components from their point-specific counterparts. The components incorporated for computing $\bar u_{{\rm V}}$ are the uncertainties in ice thickness interpolation $\bar u_{{\rm int}}$, uncertainties due to GPR measurement errors $\bar u_{{\rm gpr}}$, uncertainties in glacier area $\bar u_{{\rm area}}$ and uncertainties due to errors in surface elevation $\bar u_{{\rm surf}}$. The GPR measurement errors affect the ice volume uncertainty of glaciers both with and without GPR data, the latter via the model calibration. To distinguish between the two cases, we define $\bar u_{{\rm gpr}} = \bar u_{{\rm w/\, gpr}}$ for glaciers with GPR data, and $\bar u_{{\rm gpr}} = \bar u_{{\rm w/o\, gpr}}$ for glaciers without GPR data.

The computation of the ice volumes consists of the summation over all cells and all glaciers. When propagating the errors, we have to consider that uncertainties of individual cells are usually correlated and, for some uncertainty components, a correlation among all glaciers has to be considered as well. Details about the computation of these uncertainty components are given in Appendix D.

For propagating the individual uncertainty components to the total volume uncertainty $\bar u_{{\rm V}, i}$ of an individual glacier i, all components are again presumed to be uncorrelated, leading to an uncertainty of

(4)$$\bar u_{{\rm V}, i} = \sqrt{\bar u_{{\rm int}, i}^{2} + \bar u_{{\rm gpr}, i}^{2} + \bar u_{{\rm area}, i}^{2} + \bar u_{{\rm surf}, i}^{2} }.$$

Across the entire Swiss Alps, it is assumed that $\bar u_{{\rm int}, i}$ and $\bar u_{{\rm gpr}, i}$ are uncorrelated among different glaciers, whereas a correlation is assumed for $\bar u_{{\rm area}, i}$ and $\bar u_{{\rm surf}, i}$ (see Appendix D for further details). From the corresponding Swiss-wide values for these uncertainty components, the uncertainty of the entire ice volume in the Swiss Alps is estimated to be

(5)$$\bar u_{{\rm V}} = \sqrt{\bar u_{{\rm int}}^{2} + \bar u_{{\rm gpr}}^{2} + \bar u_{{\rm area}}^{2} + \bar u_{{\rm surf}}^{2} }.$$

Quantity $\bar u_{{\rm V}}$ refers to the uncertainty of the ice volume in entire Switzerland for the inventory year, which is 2016 (see Section 3.3). The uncertainty arising by homogenizing ice volumes to a given year y is obtained from

(6)$$\bar u_y = \sqrt{\bar u_{{\rm V}}^{2} + \left(\bar u_{{\rm hom}} \vert y-2016 \vert {A_y + A_{2016}\over 2} \right)^{2} },\; $$

where $\bar u_{{\rm V}}$ is the uncertainty computed from Eqn (5), A _y is the glacier area in the year y and A ₂₀₁₆ the one in the reference year. The homogenization uncertainty $\bar u_{{\rm hom}}$ is defined as the root-sum-square of the uncertainty of the long-term geodetic mass balances, estimated as 0.07 m a⁻¹ (Fischer and others, 2015), the effect of the uncertain extrapolation of the glacier area change on mass balance, estimated as 0.05 m a⁻¹, and the uncertainty related to the extrapolation of the year-to-year variability in mass balance to glaciers without data, estimated as 0.2 m a⁻¹.

5. Results

5.1. GPR measurements

In the framework of the current study, a total of 890 km of GPR profiles has been recorded during the years 2016 to 2019. Additionally, another 60 km were recorded in spring 2020, which have not been used for our ice thickness interpolation, but are used below in Section 5.2 for an independent validation of our ice thickness model results. Together with the GPR data of previous studies, there is now measured ice thickness along a total of around 2450 km of GPR profiles for the Swiss Alps (Table 1). These data were recorded on 251 different glaciers (red areas in Fig. 1a). Expressed in area, this is a coverage of 81% and it comprises 93% of the total glacier ice volume, computed from the volume estimate of our study. This is a remarkably high coverage. To our knowledge, such a comprehensive coverage of ice thickness observations exists for Iceland (Björnsson and Pálsson, 2020) but not for other glacierized mountain ranges at a regional scale (here almost 1000 km² of glacier area). Nevertheless, not all glaciers could be surveyed in an optimal manner, for example with regard to data quality, and with the desired density of observations. To visualize the heterogeneity of the data coverage, we computed the mean of the distances of each point on a glacier to the closest GPR point measurement, d _mean−gpr. Resulting values are listed in Table 2 and in the Supplementary material and displayed color-coded in Figure 1. The highest coverage (among the larger glaciers) exists for Glacier de la Plane Morte (A55f/03) with a distance of 57 m to the closest GPR point. The lowest coverages were found for the large glaciers Fieschergletscher (B40/07), and Gornergletscher (B56/07), the latter with the highest value of 810 m. Averaging d _mean−gpr over all glaciers with GPR data results in 187 m. The actual location of GPR survey points is displayed in gray in Figure 1b. For Fieschergletscher with its relatively low coverage, we, for example, observe that the large parts of the glacier located on steep slopes have not been surveyed.

Table 1. GPR contribution from the different sources and how many glaciers (in numbers, area and volume) are covered

^aFor sources ‘GlaThiDa’ (3.0.1) and ‘VAW’ (unpublished), see Section 2. ‘New’ refers to the data recorded in this study. Data from ‘GlaThiDa’, ‘VAW’, ‘New (2016–19)’ were used as input for the glaciological modeling, whereas data ‘New (2020)’ were recorded in spring 2020 and used for independent validation only. ‘Total’ is the sum of numbers from all sources (a specific glacier can be covered with data from different sources), and ‘Fraction’ is the percentage from the Swiss-total.

Table 2. Twenty largest glaciers by volumes, for which GPR data were used for ice thickness modeling. A complete list of glaciers with GPR data is provided in the Supplementary material

The following information is provided: areas according to the SGI2016 (A), volumes (V ₂₀₁₆), mean thickness (h _mean) and maximum thickness (h _max) of the interpolated ice thickness distributions, mean distance to closest GPR point (d _mean−gpr), total profile lengths, year of the GPR survey and the corresponding references for the GPR data. ^aGPR survey references: [1] Bauder and others (2003)^b; [2] Capt and others (2016)^a; [3] Farinotti and others (2009a) after Thyssen and Ahmad (1969); [4] Farinotti and others (2009a)^a; [5] Feiger and others (2018)^a; [6] Fischer and others (2013)^a, Huss and Fischer (2016)^a; [7] VAW (1998)^b; [8] Grab and others (2020); [9] Grab and others (2020) acquired in 2020; [10] Huss (2010)^a; [11] Huss and others (2008)^a; [12] unpublished GPR data from 2008 (ETHZ) and 2012 (UZH, UFR)^a; [13] Huybrechts and others (2008); [14] Lüthi (2000)^b; [15] Moll (2012)^a; [16] Rutishauser and others (2016)^a; [17] Sugiyama and others (2008)^a; [18] unpublished GPR data form UFR^a; [19] unpublished GPR data ETHZ^b; [20] Waechter and Roethlisberger (1982)^b; [21] Roethlisberger and Funk (1987)^b; [22] Gudmundsson (1994) and Funk and others (1994)^b; [23] Lüthi (1994)^b; [24] Sharp and others (1993); [25] VAW (2010)^b; [26] VAW (2011)^b; [27] VAW (2012)^b; [28] VAW (2014a)^b; [29] VAW (2014b)^b; [30] VAW (2017a); [31] VAW (2017b); [32] VAW (2019) ^aretrieved from GlaThiDa Consortium (2019),^bretrieved from GLAMOS (2020a).

5.2. Ice thickness distribution and glacier bed topography

The ice thickness distribution resulting from Eqn (1) is illustrated in Figure 4a for the example region of the Bernese Alps. In this region, Grosser Aletschgletscher (B36/26 in Fig. 1), the largest glacier of the European Alps, reaches a maximum ice thickness of 790 m. Subtracting the ice thickness from the DEM yields the glacier bed topography. This is shown as hillshade in Figure 4b. Parts of the glacier bed topography, which feature overdeepenings, are displayed in red colors, with the color-code indicating the overdeepening depth. For other regions of the Swiss Alps, analogous figures are provided in the Supplementary material.

Fig. 4. Ice thickness and glacier bed overdeepenings in the Bernese Alps (see box in Fig. 1): (a) ice thickness distribution. (b) Hillshade of the glacier-bed topography within SGI2016 outlines and color-coded the depth of overdeepenings. For other parts of the Swiss Alps, see the Supplementary material. Background: hillshade from the swissALTI3D after Swisstopo (2019).

5.2.1. GlaTE-, ITVEO- and MEAN performance

A comparison between modeled and measured ice thickness as well as a comparison among GlaTE, ITVEO and MEAN is shown exemplarily in Figure 5 for Kanderfirn, Tschingelfirn and Blüemlisalpgletscher (glaciers marked with IDs A55b/13, A54m/21 and A55b/02 in Fig. 1b). Profiles 3, 7 and 45 are three sections along GPR profiles (Figs 5b to d). Comparisons from all other, not shown, profiles reveal similar insights: in general, we observe a smaller amount of structural details in the modeled ice thickness distributions in comparison with the ice thickness measured with GPR (e.g. Fig. 5c). This is related to the limited resolution of the models (10 m grid size). One could argue that interpolation could be performed on an even denser grid. However, to avoid the appearance of artificial structural details within the large space in-between GPR profiles, a smoothing constraint had to be applied to both models. This counteracts the ability to fit the GPR data, independently from the grid resolution. In most cases, the models fit the measured ice thicknesses well within the measurement uncertainty, indicating that an adequate weight was given to the smoothness constraint. At the glacier margins and especially at the glacier tongues, ice thickness measured with GPR sometimes exhibit larger values than the modeled ice thickness (Fig. 5b). This is because some weight is given to the model boundary constraints, which force the models to approach zero ice thickness toward the glacier margins.

Fig. 5. Comparison of glacier bed elevations measured with GPR with modeled glacier bed elevations from MEAN, GlaTE and ITVEO. (a) Map of the region (see also Fig. 1). Location of the selected GPR profiles are shown in orange together with locations of all other profiles measured on these glaciers in black. The purple profiles are transects without GPR data. (b)–(d) Comparison of glacier bed elevations along the selected GPR-profiles. (e)–(h) Comparison of glacier bed elevations along transects without GPR data.

Intercomparing the performance of GlaTE and ITVEO, we observe only small differences between the two models within GPR cross sections. Generally, the GlaTE model converged to a slightly smoother ice thickness distributions compared to ITVEO (e.g. Fig. 5d). This is a result of the somewhat less restrictive constraints of the GlaTE model, which only targets for a fit of 95% of the cells with GPR data (Langhammer and others, 2019a), and was, thus, obtained at the expense of a somewhat poorer GPR fit. For comparing GlaTE and ITVEO at locations not covered by the GPR profiles, a transverse (A) and longitudinal (B) cross section, as well as a section across an ice divide (C), is shown in Figures 5e–g. Furthermore, a transverse cross section on Blüemlisalpgletscher (D) is shown in Figure 5h, for which no GPR data exist at all. As expected, the ice thicknesses of the two models deviate more in these cases than along GPR profiles, with differences of up to ~25% of the local ice thickness for the sections on glaciers with GPR data and up to ~35% for the glacier without any GPR data. In many cases, the ice thickness of GlaTE is larger than the one of ITVEO in the vicinity of the glacier margins, as can be seen in Figures 5f and h. This indicates a slight tendency toward more U-shaped cross sections for GlaTE and toward more V-shaped cross sections for ITVEO. For MEAN, however, pronounced features of the individual models are averaged out, yielding higher robustness against artificial structures that are not evident in the measured ice thickness.

For comparing the ice thickness distribution obtained by modeling with the measured GPR data, modeled ice thickness values were bilinearly interpolated and extracted for the exact GPR locations. Averaged over the entire Swiss Alps, the modeled values fit within the measurement uncertainty for around 85% of the GPR-locations (83% for GlaTE and 86% for ITVEO). As discussed above for the example shown in Figure 5, the main reasons for not fitting the remaining 15% of GPR data are (1) the limited resolution of the models in combination with smoothness constraints, (2) the model boundary constraints which sometimes contradict with the measured ice thickness values, and, in a few cases, (3) contradictory measurements between GPR data from different sources. The latter could not be avoided because the GPR raw data were not available to judge on the reason for the contradiction.

5.2.2. Uncertainty of interpolated ice thicknesses

For each point (x, y) of the interpolated ice thickness map, the standard uncertainty u ^±(x,  y) is estimated from Eqn (3). An example of the resulting uncertainty distribution is presented in Figure 6 for Brunnifirn, a small glacier in central Switzerland (glacier A51d/15 in Fig. 1a). This glacier was selected because some pronounced measurement-uncertainty occur in GPR profile 1. The ice thickness distribution from Brunnifirn reveals ice thicknesses up to ~140 m (Fig. 6a). The GPR data recorded along profile 1 is displayed in Figure 6b, together with its interpreted ice thickness profile shown in Figure 6c. Due to ambiguous reflections in the upper part of this profile, there are large uncertainties in the measured ice thickness. How this impacts the uncertainty of the ice thickness distribution toward smaller and larger ice thickness is presented in Figures 6d and h, respectively. The contribution from the individual error sources, interpolation uncertainty, GPR-uncertainty, and surface elevation uncertainty are displayed in Figures 6e–g and i–k. We observe the lower uncertainty bound to be strongly dominated by the GPR-uncertainty due to the large measurement uncertainty in profile 1 toward smaller ice thickness. In contrast, the GPR-uncertainty toward larger ice thickness is relatively small. Therefore, the interpolation uncertainty is the quantity dominating the upper uncertainty bound for the example of Brunnifirn.

Fig. 6. Brunnifirn as an example for estimating point-specific ice thickness uncertainties. (a) Ice thickness distribution with GPR-profile locations. (b) Example GPR-profile and (c) corresponding glacier bed interpretation with uncertainties of the GPR measurement, comparison with the ice thickness model and the model uncertainty. Total ice thickness uncertainty toward larger (d) and smaller (h) ice thickness and the corresponding contribution from interpolation- ($u_{{\rm int}}^{\pm }$), GPR- ($u_{{\rm gpr}}^{\pm }$) and surface DEM-uncertainties ($u_{{\rm surf}}^{\pm }$) in the small subpanels (e)–(g) and (i)–(k). The circles shown in (e) and (i) indicate over which distances $u_{{\rm int}}^{\pm }$, is expected to spatially correlate (see Section D).

For validating our point-specific uncertainty estimates, we compare u ^±(x,  y) with the uncertainty intervals $\Delta h_{{\rm GPR}}^{\pm }$ of measured ice thicknesses h _GPR acquired on 33 glaciers in the spring 2020 campaigns. These are data that had not been used for the ice thickness interpolation. For nine of these glaciers, older GPR data already existed and have been used for the ice thickness interpolation, whereas for the other 24 glaciers, no data were available before. We find that (i) for 39% of the points, $\Delta h_{{\rm GPR}}^{\pm }$ is entirely within u ^±(x,  y), (ii) for 25% of the points, h _GPR falls within u ^±(x,  y) but not the entire interval $\Delta h_{{\rm GPR}}^{\pm }$, (iii) for 17% of the points, h _GPR is outside u ^±(x,  y) but part of $\Delta h_{{\rm GPR}}^{\pm }$ is overlapping with u ^±(x,  y) and (iv) for 19% of the points, there is no overlap of the uncertainty ranges at all. Compliance, that is, the measurement falls within the modeling uncertainty range, is assumed for both cases (i) and (ii), whereas non-compliance is assumed for both cases (iii) and (iv). From this follows that u ^±(x,  y) accommodates a total of 39 + 25 = 64% of the measured ice thicknesses, which is in good agreement of u ^±(x,  y) being the standard uncertainty, which by definition is expected to encompass ~68% of the true ice thickness. For glaciers previously without GPR data, the compliance percentage is similar (66%) to the one for glaciers which already had GPR data (62%).

5.3. Total ice volume in the Swiss Alps

Summing over all glaciers, we find a total ice volume of 59.31 km³ for the MEAN model and the respective inventory year of all individual glaciers. The volumes resulting from the GlaTE and ITVEO models are almost identical, being 59.30 and 59.32 km³, respectively. This high agreement is accidental. In fact, the two models differ slightly more when differentiating between glacier size classes. After homogenization, we find a total ice volume in the Swiss Alps for the year 2016 of V ₂₀₁₆ = 58.8 ± 2.5 km³. For the magnitudes of the uncertainty components for the total ice volume we find $\bar u_{{\rm int}}$ = 0.24 km³, $\bar u_{{\rm w/\, gpr}}$ = 0.14 km³, $\bar u_{{\rm w/o\, gpr}}$ = 2.44 km³, $\bar u_{{\rm area}}$ = 0.30 km³ and $\bar u_{{\rm surf}}$ = 0.54 km³. Substituting into Eqn (5), it results an uncertainty of $\bar u_{{\rm V}}$ = ±2.54 km³ or of 4.3%, respectively.

To investigate, how the uncertainty of glaciers of different size and with/without GPR data contribute to the total ice volume uncertainty (for individual glaciers see Fig. S13 in the Supplementary material), the cumulative uncertainty $\bar u_{{\rm V}}$ is consecutively calculated from Eqn (5) for all glaciers with average ice thickness up to a certain size (Fig. 7a). The cumulative ice volume for glaciers with and without GPR data, and the cumulative glacier area, are displayed in Figure 7b. It becomes evident that the predominant part of the volume uncertainty stems from glaciers with average ice thicknesses <50 m, although they only account for around one quarter of the total ice volume. The main factor controlling the uncertainty of these smaller glaciers is $\bar u_{{\rm w/o\, gpr}}$, that is the uncertainty of calibrating the ice thickness models for glaciers without GPR data. For all other glaciers, for which GPR data have been acquired, $\bar u_{{\rm surf}}$ is the most important contribution to the uncertainty, followed by $\bar u_{{\rm area}}$, $\bar u_{{\rm int}}$ and $\bar u_{{\rm w/\, gpr}}$.

Fig. 7. (a) Overall ice volume uncertainties and share of the different uncertainty components computed by cumulating all glaciers according to their average ice thickness. (b) Cumulative ice volume for glaciers with and without GPR data, and cumulative glacier area.

6. Discussion

6.1. Ice volumes from earlier studies

The total ice volume in the Swiss Alps resulting from our study is V ₂₀₁₆ = 58.7 ± 2.5 km³ and, when extrapolating the ice volume to 2020, V ₂₀₂₀ = 52.9 ± 2.7 km³. Earlier estimates of the Swiss glacier volume were presented in various studies based on the glacier inventories of 1973 (Muller and others, 1976; Maisch, 2000), 1999 (Farinotti and others, 2009a; Linsbauer and others, 2012) and 2010 (Fischer and others, 2015). Furthermore, Farinotti and others (2019) estimated the total ice volumes in the European Alps for the year 2003 and Haeberli and others (2019) for the year 2018. The corresponding volume estimates are displayed in Figure 8. Based on empirical relationships between glacier area and mean ice thickness, Muller and others (1976) estimated a volume of 67 km³ and Maisch (2000) of 74 km³ for the year 1973. Estimations based on measured ice thickness for a subset of glaciers combined with modeled ice thickness have been obtained by Farinotti and others (2009a) with 75 ± 9 km³ and Linsbauer and others (2012) with 65 ± 20 km³ both for the year 1999, and by Fischer and others (2015) with 60 km³ referring to the year 2010.

Fig. 8. Temporal extrapolation of the total ice volume in Swiss glaciers for the period 1973–2020 (black curve with gray uncertainty range from Eqn 6) based on the overall volume in 2016 inferred in the current study (red asterisk), see Section 3.3 for details. The time series of overall glacier volume allows direct comparison to previous ice volume estimates for the Swiss Alps (green).

First of all, it is noteworthy that our estimated uncertainty is clearly smaller than the one by Farinotti and others (2009a) and Linsbauer and others (2012), which we attribute to the larger amount of ice thickness data used compared to previous studies. Furthermore, when homogenizing our results to the corresponding years of the previous estimates we find ice volumes of 65.2 ± 2.8 km³ (2010), of 77.2 ± 4.6 km³ (1999) and of 94.0 ± 10.9 km³ (1973) (Fig. 8). This indicates that the early studies by Muller and others (1976) and Maisch (2000) likely underestimated the volumes that were present in the 1970s in comparison with our results. High agreement was found with the estimate by Farinotti and others (2009a) with a discrepancy of only ~2 km³. Also for the study by Linsbauer and others (2012) the error bars overlap, although their ice volume estimate is substantially lower by ~12 km³. The estimate by Fischer and others (2015) for the year 2010 is close to our estimate but again slightly lower by ~5 km³, which might be due to the somewhat underestimated glacier area in their inventory (Fischer and others, 2014) in comparison with the SGI2016 (GLAMOS, 2020b). Huss and Farinotti (2012) estimated that Swiss glaciers contribute with 58.7% to the total ice volume of the European Alps. Based on this percentage, the volume estimate by Farinotti and others (2019) is equivalent to 75 ± 18 km³ for the year 2003, which closely fits our estimate for the corresponding year. For Haeberli and others (2019) it is equivalent to 49 ± 15 km³, which is lower by ~6 km³ in comparison with our results.

It is apparent from Figure 8 that from the newer studies, the ones by Linsbauer and others (2012) and by Haeberli and others (2019) tend to underestimate the ice volumes with respect to other recent estimates and the result of our assessment. It is difficult to investigate whether such differences are related to the ice thickness observations that were available for constraining the individual estimates, or whether they can be attributed to the different type of ice thickness model used. Model types differing from GlaTE and ITVEO are for instance models based on the shallow ice approximation (e.g. Linsbauer and others, 2012) or 2-D approaches based on the continuity equation (e.g. Morlighem and others, 2011) and on the parallel ice sheet model (e.g. Van Pelt and others, 2013). Although ITMIX2 (Farinotti and others, 2021) does not provide a direct answer in our case, no significant differences in performance or biases regarding computed ice thickness between models of different types were found. Overall, ITMIX2 revealed no drift in the model-derived ice thicknesses when restricting the input to fewer observational data. This was the case for both, individual models and the model ensemble. The absence of such a drift indicates that our ice volume estimate at the mountain-range scale is relatively robust with regard to the model types. Therefore, we expect only minor discrepancies in the total ice volume if another model type would have been constrained with the comprehensive set of ice thickness data available today. It is probable, however, that the model type has a larger impact on the volume estimate of individual glaciers or the ice thickness distribution at the local scale. Within the ITMIX2 experiment, GlaTE and ITVEO were part of the group of models giving preference to the reproduction of observational data over internal model consistency. Therefore, including a model that emphasizes internal model consistency might improve our results for glaciers with sub-optimally distributed point observations. Such a model is for instance the one presented by Van Pelt and others (2013), which showed a strong performance in ITMIX2.

As introduced in Section 3.3 and presented in Figure 8, our homogenization approach also allows extrapolating the overall ice volume of Swiss glaciers backward in time. Our results indicate a volume of 65.2 ± 2.8 km³ in 2010, 76.5 ± 4.4 km³ in 2000 and 95.1 ± 9.0 km³ in 1980. This indicates that ~0.9 km³a⁻¹ ice were lost between 1980 and 2000, 1.1 km³a⁻¹ between 2000 and 2010 and 1.2 km³a⁻¹ between 2010 and 2020. With reference to the earlier ice volumes, it is ~43% of the Swiss glacier volume that has been lost over the last 40 years, 29% within the last 20 years and around 17% within the last 10 years. These relative ice volume changes confirm the recent acceleration of glacier mass loss presented for example by Zemp and others (2015) and Beniston and others (2018).

6.2. GPR data coverage

The measured ice thickness forming the basis of our study will be publicly available. For some of the older data this is already the case through database GlaThiDa (GlaThiDa Consortium, 2019) and for all the new data and some data from older studies it is made available in the framework of the study presented here. As displayed in Figure 1, these data cover the overwhelming part of the glacierized areas in the Swiss Alps. Our uncertainty analysis, however, shows that the few glacier areas remaining without GPR data contribute the most to the uncertainty in the Swiss-wide total ice volume due to the uncertainty in model calibration for those glaciers. Therefore, performing GPR-campaigns on the largest glaciers remaining without GPR data to date, would be the most effective option for further reducing this uncertainty. The largest glacier (in volume) without GPR data is Dammagletscher (marked with ID A51f/10 in Fig. 1) with an area of 3.8 km² and an estimated volume of 0.12 km³. A list with the largest glaciers without GPR data is given in the Supplementary material.

The uncertainty arising from the interpolation procedure between the sparse grid of GPR-measurement locations contributes relatively little to the overall uncertainty of our total ice volume estimate. Therefore, we conclude that no further refinement of the GPR coverage would be required for future volume estimates or similar studies focusing on scales beyond the one of individual glaciers. Recommendations about survey layouts at such scales was deduced from our earlier study (Langhammer and others, 2019a) where we proposed cost-optimized GPR survey layouts.

The interpolation uncertainty contributes significantly to the point-specific uncertainty in ice thickness as, for example, shown in Figure 6d. For specific studies on smaller scales, it might therefore be recommendable to acquire GPR data on denser grids than we have typically done. An example for such a survey is a detailed study that we have performed on the tongue of Oberaletschgletscher, for which data were acquired on a dense grid with around 70 m profile-interspacing as shown in Figure 1 (Oberaletschgletscher marked with ID B36/01).

6.3. Ice thickness distribution and glacier bed topography

The ice thickness maps and the glacier bed topography maps resulting from our study are made publicly available as well. We anticipate these data to be of value for various further glaciological investigations or studies in related fields. For example, the study of Gabbi and others (2012) investigated river runoff projections using a glacio-hydrological model that requires detailed information about the ice thickness distribution as a function of elevation and slope exposition. It was frequently stressed that more accurate ice thickness data are the most important factor for reducing the uncertainty in long-term projections of glacier retreat and corresponding changes in the runoff regime (e.g. Huss and others, 2014), for example to obtain the precise timing of maximum glacier melt volumes (Schaefli and others, 2019), or to better predict future hydropower revenue (Gaudard and others, 2016). From our data, for example, the total ice volume for individual river catchments and across different elevation bands can be calculated (see Supplementary material). It is important to note here that we homogenized the ice volumes but not the ice thickness distributions. The latter are given for the inventory years defined by the DEM (swissALTI3D and r2019) and thus exhibit some heterogeneity. For this reason, we also deliver the meta data, which provides the inventory year in the form of raster data (Supplementary Fig. S11).

The maps of glacier bed topography provide information on the future valley morphology forming after glacier retreat, which is expected to further accelerate in the coming decades (Zekollari and others, 2019). Such a map, for example, contains information on glacier bed overdeepenings, which can be further analyzed with regard to their potential for the formation of future proglacial lakes. Based on the glacier bed topography and the distribution of today's ice volume, various studies have aimed at predicting the size and the formation time of such lakes in the Alps (Frey and others, 2010; Linsbauer and others, 2012; Gharehchahi and others, 2020; Magnin and others, 2020; Shugar and others, 2020). Such information is important since future lakes might change the hazard situation in the corresponding regions, or can be large enough to be of economic interest for tourism, fresh water supply or in the context of hydropower (Haeberli and others, 2016). We have thus analyzed the glacier bed map resulting from our study with regard to locations and volumes of overdeepenings. For the example region of the Bernese Alps, the location of anticipated overdeepenings is displayed in Figure 4b with the color code indicating their depth. In comparison with the earlier study by Linsbauer and others (2012), there is a good agreement with regard to the locations where the most prominent overdeepenings occur. Computing the total volume of these overdeepenings across entire Switzerland results in a volume of around 1 km³ and delivers an estimate of future lake volumes. More than half of this lake volumes occur in the glacier beds of only four glaciers (Supplementary Fig. S9c). For the same glaciers except for Rhonegletscher, also Linsbauer and others (2012) observed large overdeepening volumes. The total volume, however, is substantially smaller than the lake volume of ~2 km³ estimated by Linsbauer and others (2012). The discrepancy can to some degree be explained by the different modeling technique used by Linsbauer and others (2012) in comparison with the one used here. Nevertheless, since our results are based on substantially larger amounts of measured ice thickness data, it is likely that the total volume of lakes forming in the future in the Swiss Alps is smaller than previously assumed.

When referring to our glacier bed topography, it is important to note that it might also consist of unconsolidated sediment rather than bedrock. Reflections in the GPR images interpreted as the glacier bed, typically originate from the uppermost horizons with a pronounced GPR impedance contrast. In case of porous basal sediment layers saturated with liquid water or of a liquid water film on top of a well consolidated sediment deposit, such pronounced impedance contrasts occur at the base of the ice rather than at the base of the sediment layers. Thus, we expect that our glacier bed topography in general represents the topography at the base of the ice.

7. Conclusions

We have measured the ice thickness of all larger and most medium-sized (≥0.02 km³) glaciers of the Swiss Alps that had not been surveyed before, and we have refined the measurement grids on some of the glaciers with only partial surveys from previous studies. Together with the pre-existing data, measured ice thickness data are now available for 251 glaciers, which account for 81% of the glacierized area in the Swiss Alps. Altogether, these glaciers contain 93% of the total ice volume. With this comprehensive data coverage, our dataset has a great potential for various further investigations.

The measured ice thickness data, available only on the sparse grid of GPR profiles, were interpolated using two independent glaciological modeling techniques, GlaTE and ITVEO, together with the glacier outlines from the Swiss Glacier Inventory 2016 and the digital elevation model swissALTI3D-r2019. The use of more than one algorithm leads to a higher robustness of our results, in particular against interpolation artifacts. We obtained continuous ice thickness and glacier bed topography maps for all glaciers in the Swiss Alps at 10 m resolution, and estimated a total ice volume of 58.7 ± 2.5 km³ for the year 2016 and of 52.9 ± 2.7 km³ for 2020. Thanks to the much more complete coverage of measurements, we were able to estimate the total ice volume, the ice thickness distribution and glacier bed topography with substantially lower uncertainty than in previous studies. The overall uncertainty is mainly controlled by the uncertainty in the ice volume of glaciers remaining without GPR data, which accounts for only 7% of the total ice volume. The volumes of these glaciers were estimated by extrapolating the model calibration from glaciers with GPR data to glaciers without GPR data.

Our results are accessible as the SwissGlacierThickness-R2020 data collection (Grab and others, 2020) for further use. This data collection also includes the ice thickness point data, which will be additionally available through a future release of GlaThiDa (Welty and others, 2020). We anticipate these data to be of value for various further studies, as, for example, investigations of river runoff projections, or for predicting locations and volumes of proglacial lakes that will occur in the course of glacier retreat. An analysis of our glacier bed topography maps, for example, revealed that the volume of such lakes forming in the future in the Swiss Alps is probably smaller than it was assumed from earlier studies.