Introduction
Mountain glaciers and ice caps (GICs) are important contributors to current sea-level rise (SLR) (Reference SolomonSolomon and others, 2007) and are likely to continue to be so during the 21st century (Reference Radic and HockRadić and Hock, 2011), in spite of the growing contribution from the ice sheets of Greenland and Antarctica (Reference Rignot, Velicogna, Van den Broeke, Monaghan and LenaertsRignot and others, 2011). For how long, and to what extent, GICs will contribute to SLR depends on their total ice volume, which is not accurately known because of the scarcity of direct observations of glacier ice thickness. Because of this, the global estimates of the volume of GICs are based on either volume/area scaling relationships (Reference Radic Vand HockRadić and Hock, 2010; Reference GrinstedGrinsted, 2013) or on simple physically based approaches relating thickness distribution to glacier geometry and dynamics (Reference Huss Mand FarinottiHuss and Farinotti, 2012). Their estimates, however, vary considerably.
Svalbard, situated in the Atlantic sector of the Arctic (76–818 N, 10–338 E), is highly vulnerable to climate change (Reference Hagen, Melvold, Pinglot and DowdeswellHagen and others, 2003). Even though Svalbard is among the best-studied sectors of the Arctic, ice thickness measurements are only available for selected glaciers. The estimates of the total ice volume of Svalbard glaciers are therefore based on the scaling or physically based approaches mentioned above, or on empirical relationships specific to Svalbard glaciers, derived from the scarce and limited-accuracy ice thickness measurements that were available in the early 1980s (Reference YuYa and ZhuravlevMacheret and Zhuravlev, 1982) and the early 1990s (Reference Hagen, Liestøl, Roland, J and ørgensenHagen and others, 1993).
The results presented in this paper are part of ongoing work within the ‘Sensitivity of Svalbard glaciers to climate change’ (SvalGlac) project. One of its major objectives is to obtain a reliable estimate of the total ice volume of Svalbard glaciers by deriving a new volume/area relationship specific to Svalbard glaciers, based on accurate ice volume estimates from ground-penetrating radar (GPR) data. Here we use GPR data collected during several field campaigns within the period 1999–2012 to calculate the ice volume of ten glaciers in western Norkdenskiöld Land, central Spitsbergen, Svalbard (Fig. 1), together with associated error estimates. This estimation involves compound errors from different sources, which often are not properly quantified in the literature. In addition to the volume estimates, this paper provides an insight into how to deal with error estimates when computing glacier volume from GPR-retrieved ice thickness data.
Fig. 1. Location of the studied glaciers and layout of GPR profiles.
Field Data
The GPR data used in this study correspond to ten glaciers in Nordenskiöld Land, collected during several field campaigns between 1999 and 2012, all of them carried out during early spring, before the onset of strong surface melting. The different types of radar equipment, with their central frequency and the total length of radio-echo sounding (RES) profiles, are summarized in Table 1. The layout of GPR profiles is displayed in Figure 1 , which also includes the 2007 glacier boundaries (Arendt and others, 2012; König and others, in press), modified using Landsat images to match the area at the time when the glaciers were echo sounded. For the profiling, transmitting and receiving antennas were arranged coaxially along the profiling direction (parallel end-fire), to minimize direct coupling between antennas and also the reflections from the glacier side-walls, as most of the profiles were transverse to the glacier centre line (Reference Navarro, Eisen, Pellikka and ReeseNavarro and Eisen, 2010).
Table 1. GPR surveys. Ramac/GPR is manufactured by Mala Geoscience. VIRL2, VIRL6 and VIRL7 are GPR systems made in-house (described by Reference Vasilenko, Sokolov, YuYa, Glazovsky, Cuadrado and NavarroVasilenko and others, 2002; Reference Berikashvili, Vasilenko, Macheret and SokolovBerikashvili and others, 2006, and Reference Vasilenko, Machío, Lapazaran, Navarro and FrolovskiyVasilenko and others, 2011, respectively)
Methods
GPR data processing
The radar data were processed using the commercial software packages RadExPro, by GDS Production (Reference Kulnitsky, Gofman and TokarevKulnitsky and others, 2000), and ReflexW (http://www.sandmeier-geo.de/), and the main processing steps consisted of bandpass filtering, amplitude correction, deconvolution and migration. For the time-to-depth conversion we used a constant radio-wave velocity (RWV) of 170m μs_1. This RWV value was chosen on the basis of our own measurements by the common-midpoint (CMP) method (Reference YuYa, Moskalevsky and VasilenkoMacheret and others, 1993) on glaciers in this region (Reference Navarro, Glazovsky, YuYa, Vasilenko, Corcuera and CuadradoNavarro and others, 2005; Reference Vasilenko, Glazovsky, Lavrentiev, YuYa and NavarroVasilenko and others, 2006) or in neighbouring regions in Svalbard (Reference JaniaJania and others, 2005), carried out during early spring. Such a high RWV is typical of measurements made before the onset of strong melting; measurements made during warmer periods have substantially lower RWVs, due to the increased water content in temperate ice. A RWV of 170 m μs- 1 has also been used in other GPR studies in Svalbard (Reference Bælum and BennBælum and Benn, 2011; Reference SaintenoySaintenoy and others, 2013).
Data consistency
To properly interpolate a continuous surface from point data a coherent (self-consistent) source dataset is required. Crossover analysis is used to assess both internal accuracy and consistency between datasets (Reference Bamber, Layberry and GogineniBamber and others, 2001). Since all data contain errors, the ice thickness values measured at any crossover point will, in general, be different for the two intersecting radar lines. The value and spatial distribution of the crossover errors provide insight into the magnitude and source of the errors in the data (Reference Retzlaff, Lord and BentleyRetzlaff and others, 1993). Differences at the crossover points can be caused by inaccurate positioning of the data, picking errors or lack of (or improper) migration of any of the intervening profiles. In our analysis, the highest crossover differences appeared for crossovers involving radar profiles parallel and close to steep valley side-walls. In such cases, the radar often receives reflections from (and perpendicular to) the slopes of the valley walls rather than from the point underneath the radar position. The standard migration procedure only corrects the ice thickness for bed slope along the profiling direction. Consequently, it works properly for profiles perpendicular to the side-walls, but not for profiles parallel to them, and hence a large difference often appears at such crossovers. To avoid these problems, we discarded profiles parallel and close to the side-walls.
The ice thickness maps were constructed assuming that: (1) the ice thickness is zero at the lateral margins of the glacier, at the glacier front of land-terminating glaciers and at the contact points with nunataks; (2) if no thickness measurements are available at the other side of an ice divide, the ice thickness at the divide is extrapolated by a second-degree polynomial passing through three selected measurement points close to the divide and roughly aligned in the longitudinal (glacier-flow) direction; (3) the ice thickness near the valley side-walls decreases parabolically to zero at the side-wall along the direction perpendicular to it; and (4) the ice thickness near the front of land-terminating glaciers decreases linearly to zero at the front, a reasonable assumption for receding land-terminating glaciers (Reference Cuffey and PatersonCuffey and Paterson, 2010), such as those in this study (with the exception of Fridtjovbreen).
Interpolation
Ice thickness measurements were interpolated into a regular grid using the Gridfit routine (Reference D’ErricoD’Errico, 2006), available for MATLAB. This method has the capability to work with sparse or irregular data, and finds a smooth surface that is consistent with the data. This approach is better than implementing smoothing after interpolation, where smoothing is disconnected from the data. Gridfit uses an approximation method closely related to thin plate splines approximation, as opposed to interpolation, so that it is robust against noise and outliers (Reference Bohorquez and DarbyBohorquez and Darby, 2008). The Gridfit routine has different user-selectable options, such as the interpolation method, the smoothing factor or the regularization method. We used the triangular interpolation method, leaving the default values for the remaining parameters. Further details about the algorithms used by Gridfit and its efficiency can be found at MATLAB Central File Exchange (Reference D’ErricoD’Errico, 2006).
Volume calculation in brief
The procedure we followed to derive the glacier volumes consisted of the following steps. (1) Collecting a coherent GPR dataset, as discussed above. (2) Processing of GPR data following the methodology described above. (3) Calculating ice depths along the profile at ∼1 . 5m intervals. (4) Interpolating the data into a regular grid using MATLAB. (5) Adjusting the 2007 glacier boundaries (Arendt and others, 2012; König and others, in press) to the year of the RES using Landsat images (this and the following step were performed using ArcGIS v9.3). (6) Clipping the resulting regular grid to the glacier boundary to build the final ice thickness map for each glacier. (7) Calculating the glacier volume, summing up every gridcell thickness value multiplied by the gridcell area.
Error in volume estimate
We attempt to estimate the error in volume, using error propagation, from the separate error estimates for its main components. Errors in ice thickness and errors in area can be considered independent of each other, so their combined effect is approximated by root-mean-square (rms) summation
Measurement error
The measurement error encompasses errors inherent to the ice thickness measurements using GPR (GPR errors) and errors due to inaccurate GPS positioning (GPS errors). We assume that all ice thickness data are affected by random errors. GPR errors encompass both instrumental errors and data-processing errors. The instrumental errors include the range resolution (which depends on the radar frequency) and the data-digitalization error (which depends on the sampling frequency). The range resolution is not an error in itself, but acts as a lower bound for the thickness measurement error. The data-processing errors include, among others, those related to lack of (or improper) migration, improper bed reflection picking or the errors incurred by the assumptions in the time-to-depth conversion (e.g. constant RWV). We assume that picking errors are negligible compared to the other errors. To estimate the GPR measurement error we followed the technique described by Reference Navarro, Eisen, Pellikka and ReeseNavarro and Eisen (2010), which takes the error in two-way travel time as the inverse of the central frequency of the radar system, and assumes a typical relative error in RWV of ∼2%.
Stand-alone GPS systems were used for positioning. For the field surveys carried out after 2000, a typical accuracy of 5 m in horizontal positioning can be assumed (Reference Sharma and BanerjeeSharma and Banerjee, 2009). However, for the 1999 echo sounding on Aldegondabreen the accuracy of the GPS horizontal positioning is expected to be much lower (errors of up to 100 m) (Reference Parkinson, Spilker, Axelrad and EngeParkinson and others, 1996). To estimate how much thickness change could be expected from inaccuracies in the horizontal positioning, we built a small-scale experimental variogram, assessing the thickness differences between measurement points separated by a distance equal to the GPS precision (100 m for Aldegondabreen, 5 m for the other glaciers). This quantifies the thickness error incurred by assigning to a given point the ice thickness corresponding to a different position, due to inaccurate GPS horizontal positioning. This difference will be larger in zones of steep bedrock.
Interpolation error
GPR data are usually anisotropically distributed over the glacier surface, densely sampled along the radar tracks but widely spaced between distinct tracks. Well-known methods, such as cross-validation or jack-knifing techniques (Reference DavisDavis, 1987; Reference Isaaks and SrivastavaIsaaks and Srivastava, 1989), have often been used to quantify the accuracy of surface interpolation algorithms. However, these errors are not good estimators of the interpolation error when dealing with GPR-like data, because of their sparse and uneven distribution. Errors will not be representative for the glacier areas not covered by GPR profiles. Our method calculates an average interpolation error following the rationale of cross-validation techniques, but taking into account the variance of the error with the distance to the nearest neighbour. We first construct a function relating the interpolation error at any given point to the distance to the nearest GPR-measured data point, and then calculate a glacier-wide interpolation error by averaging the interpolation errors computed at each gridpoint.
Error in area
The error in area has two main sources: the uncertainty in identifying the boundary of the glacier (which we will refer to as uncertainty in boundary delineation) and, once a given boundary has been assumed, the error incurred in fitting the glacier boundary to the regular grid (which we will denote pixellation error). The uncertainty in boundary delineation arises from the overestimation of the area in the accumulation zone, due to snow-covered terrain, and the underestimation incurred in the ablation zone when glacier ice near the margins is covered by debris. In our study area, debris cover is fairly common, taking into account that most of the glaciers are small land-terminating mountain glaciers with steep valley walls. To account for the uncertainty in boundary delineation we have assumed an error in area of 8%, typical of glacier inventories for this region (personal communication from C. Nuth, 2013). The pixellation error has a negligible impact on the volume computation compared with the boundary delineation error assumed in this study.
Results and Discussion
The computed glacier areas and volumes are shown in Table 2, together with the detail of the estimations for the different sources of the error in volume. The ice thickness maps are displayed in Figure 2. The individual areas, volumes and average ice thickness lie within 2.6–50.4 km2, 0.08–5.54 km3 and 29–108 m, respectively. The maximum recorded ice thickness, 265 ± 16 m, corresponds to Fridtjovbreen, which also has the largest average thickness (108 ± 1 m). The total area and volume of the ensemble are 116:06 ±4:53 km2 and 10:439 ± 0:373 km3, respectively.
Table 2. Estimated area based on the boundaries of the Randolph Glacier Inventory (RGI), volume, average thickness and maximum thickness for each glacier, with their estimated errors, and main components of the error in ice thickness
Fig. 2. Ice thickness maps of the studied glaciers.
The relative errors in volume are generally <10%, except for Aldegondabreen (due to the large GPS error) and Vestre Dahlfonna (due to the large interpolation error compared with its mean thickness). These errors are larger than those usually found in the literature, which we believe are unrealistically small because they are often based on the assumption that all errors are independent.
We have taken RWV to be 170(±2%) m μs– 1 for the time-to-depth conversion, on the basis of field measurements in the region, with the 2% accounting for random errors, because the RWV is known to show both space and time variations, due to changes in material density or water content (Reference Navarro, Eisen, Pellikka and ReeseNavarro and Eisen, 2010). Both the RWV measurements and the GPR profilings reported in this paper were carried out during early spring, before the onset of surface melting, which minimizes the magnitude of the temporal variations of RWV. Regarding the spatial variations, our analysis shows that, for the glaciers studied and the early springtime, the column-averaged RWVs in the accumulation and ablation zones do not differ substantially. We attribute this to the fact that, despite the higher RWV in firn (∼190 m μs–1), the firn layer is, in general, not too thick and it is underlain by temperate ice with typical velocities of 165 m μs– 1 or lower, so the column-averaged velocity in the accumulation zone is close to the ∼170 m μs– 1 measured at some CMP locations on the upper ablation zone (Reference JaniaJania and others, 2005; Reference Navarro, Glazovsky, YuYa, Vasilenko, Corcuera and CuadradoNavarro and others, 2005). In addition to random errors, systematic errors could also intervene. The most likely one would be the bias introduced by selecting a constant RWV systematically above (or below) the real values almost everywhere in the glacier. The relative error in volume introduced by a such a hypothetical bias in RWV can be readily estimated: it is equal to the relative error in RWV, because ice thickness is calculated as half the product of two-way travel time and RWV, and volume is computed as area multiplied by thickness.
It is of interest to analyse how well the global and regional volume/area scaling relationships perform for this region. Table 3 summarizes the results of the comparison between the measured volumes and those computed from different scaling approaches. As expected, all relationships perform better when applied to the entire set of glaciers rather than to individual glaciers. The largest global discrepancies are those by Reference Hagen, Liestøl, Roland, J and ørgensenHagen and others (1993), in spite of being specific for Svalbard, and by Reference Radic Vand HockRadić and Hock (2010), which has been suggested to overestimate the volumes (Reference GrinstedGrinsted, 2013). At an individual level, the volume overestimate is largest for the thinnest glaciers (Gleditchfonna, Tungebreen and Vestre Dahlfonna). The relationships of Reference Chen and OhmuraChen and Ohmura (1990) and Reference Bahr, Meier and PeckhamBahr and others (1997) give a good overall result and also behave reasonably well at the individual level, while others producing good results globally show large discrepancies at the individual level (e.g. Reference YuYa and ZhuravlevMacheret and Zhuravlev, 1982). The fact that the Reference YuYa and ZhuravlevMacheret and Zhuravlev (1982) relationship is based on a sample of Svalbard ice masses which includes some glaciers analysed here could explain its good behaviour when applied to our dataset. These results should be treated cautiously, as one of the glaciers, Fridtjovbreen, has a volume nearly as large as the sum of all the others, implying a strong bias. In fact, if this glacier is removed from the dataset, the results change substantially, although Reference Hagen, Liestøl, Roland, J and ørgensenHagen and others (1993) and Reference Radic Vand HockRadic and Hock (2010) still produce the largest overestimates. Considering both cases (Fridtjovbreen included in and excluded from the dataset), the scaling law by Reference GrinstedGrinsted (2013) using the World Glacier Inventory (WGI)/ Global Land Ice Measurements fom Space (GLIMS) inventory gives the best overall result. We note that Fridtjovbreen, in spite of its surge in the 1990s (Reference Murray, Luckman, Strozzi and NuttallMurray and others, 2003), does not show an anomalous deviation from the values provided by the area/volume relationships.
Table 3. Comparison of GPR-retrieved volumes, V GPR (km3), with volumes computed using area/volume scaling relationships, V sc, given as per cent differences (100 x (V sc – V GPR)=V GPR)
Of the glaciers analysed in this paper, volume estimates from ice thickness data retrieved from RES with a dense coverage of the entire glacier surface have only been reported for Aldegondabreen in 1999 (Reference Navarro, Glazovsky, YuYa, Vasilenko, Corcuera and CuadradoNavarro and others, 2005). The volume reported was 0.558 km3, without an accompanying error estimate. Our present volume estimate, based upon the same GPR field data, is ∼16% smaller. This discrepancy is mostly introduced by the removal, in our volume estimates, of the radar profiles parallel and close to the very steep side-walls of Aldegondabreen. This was recommended by the results of our crossover analysis. Removing such profiles resulted in smaller glacier thickness near both the northern and southern side-walls. In these areas, a large amount of seasonal snow remains accumulated in early spring and, given the high RWV of dry snow (typically 215-255 m μs– 1; e.g. Reference Harper and BradfordHarper and Bradford, 2003), the glacier thicknesses calculated for these zones by Reference Navarro, Glazovsky, YuYa, Vasilenko, Corcuera and CuadradoNavarro and others (2005) are expected to be overestimated, leading to a larger total volume. Volume data reported previously for the remaining glaciers analysed in this paper were based either on volume/average thickness relationships (Reference Hagen, Liestøl, Roland, J and ørgensenHagen and others, 1993) or on estimates using airborne RES data along the glacier centre lines and assuming a parabolic shape for the glaciers (Reference Macheret and ZhuravlevMacheret and Zhuravlev, 1980, Reference YuYa and Zhuravlev1982; Macheret, 1981).
The presence or absence of scattering in the radargrams provides a way to interpret the glacier ice as temperate or cold (Reference Navarro, Eisen, Pellikka and ReeseNavarro and Eisen, 2010), allowing inference of the hydrothermal structure of the main glaciers of western Nordenskiold Land. The hydrothermal structure of Fridtjovbreen, prior to its surge in the 1990s (Reference Murray, Luckman, Strozzi and NuttallMurray and others, 2003), was interpreted from airborne (Reference Macheret and ZhuravlevMacheret and Zhuravlev, 1980; Reference Dowdeswell, Drewry, Liestøl and OrheimDowdeswell and others, 1984a) and ground-based (Reference Glazovskiy, Macheret, Moskalevskiy and JaniaGlazovskiy and others, 1991) RES data. However, neither the Russian (Reference Macheret and ZhuravlevMacheret and Zhuravlev, 1980, Reference YuYa and Zhuravlev1982; Macheret, 1981) nor the British (Reference Dowdeswell, Drewry, Liestøl and OrheimDowdeswell and others, 1984a,Reference Dowdeswell, Drewry, Liest, øl and Orheimb) echo soundings revealed the poly-thermal structure of Austre Grønfjordbreen shown in our GPR data. Neither was the Russian airborne echo sounding able to detect the polythermal structure of Aldegondabreen, though it was later inferred from ground-based GPR (Reference Navarro, Glazovsky, YuYa, Vasilenko, Corcuera and CuadradoNavarro and others, 2005). GPR data from the campaigns from 2010 onwards have revealed the polythermal structure of Tavlebreen (Reference Lavrientiev, YuYa, Holmlund and GlazovskyLavrientiev and others, 2010), Austre and Vestre Grønfjordbreen and Tungebreen (Reference &&&&Lavrientiev, YuYa and GlazovskyLavrientiev and others, 2011), as well as the cold structure of Gleditchfonna (Reference &&&&Lavrientiev, YuYa and GlazovskyLavrientiev and others, 2011). GPR data reveal that Erdmanbreen, Austre Dahlfonna and Vestre Dahlfonna are also polythermal.
Conclusions
The main conclusions of our analysis are:
The total area and volume of the ensemble of ten western Nordenskiold Land glaciers considered in this study are 116:06 ± 4:53 km2 and 10:439 ± 0:373 km3, respectively, while the individual areas, volumes and average ice thickness lie within 2.6-50.4 km2, 0.08-5.54 km3 and 29-108 m, respectively. The maximum ice thickness (265 ± 16 m) was recorded on Fridtjovbreen.
The total volume of these glaciers calculated using volume/area scaling relationships is overestimated by up to 35% compared to our results. The scaling approach results are biased by the large volume of Fridtjovbreen, which accounts for 45% of the total volume. In general, glaciers with the smallest average ice thickness provide the worst scaling results, with some volumes overestimated by >100%.
On the basis of the pattern of scattering in the radargrams, we suggest that Aldegondabreen, Vestre and Austre Dahlfonna, Erdmanbreen, Fridtjovbreen, Vestre and Austre Grønfjordbreen, Tavlebreen and Tungebreen are polythermal glaciers, while Gleditchfonna seems entirely cold.
Acknowledgements
This research was supported by grant EUI2009-04096 from the Spanish EuroResearch Programme (a part of the SvalGlac Project, PolarCLIMATE Programme of the European Science Foundation), grants CTM2008-05878/ANT and CTM2011-28980 from the Spanish National Plan for R&D, and grants 10-05-00133-a and 11-05-00728-a from the Russian Fund of Basic Research. The comments by the Scientific Editor, John Woodward, and an anonymous referee contributed to improve the manuscript.
