GPR DATA PROCESSING AND VOLUME CALCULATION

The radar data were processed using the commercial software packages RadExPro, by GDS Production (Kulnitsky et al., 2000), and ReflexW (Sandmeier Scientific Software, 2012). The main processing steps consisted of bandpass filtering, amplitude correction, deconvolution, and migration. For Ramac/GPR data, predictive deconvolution was used (e.g., Yilmaz, 2001), whereas for VIRL GPR data, deconvolution was not used because our pulse duration is small (∼25 ns) and thus there is no need to shorten it. The migration algorithm applied was Stolt F-K (e.g., Yilmaz, 2001). Its suitability was verified by the proper collapse of the diffraction hyperbolae. For improving the detection of the zero time of the transmitted pulse and of the bed reflections, a Hilbert transform (e.g., Taner et al., 1979) was applied upon filtering (which results in small picking errors) of the order of the sampling period (usually within 2.5–5 ns, depending on the radar campaign). For the time-to-depth conversion we used a constant radio-wave velocity of 166 m µs for the largest glaciers (Hansbreen, Paierlbreen, Recherchebreen, and Austre Torellbreen), whereas 168 m µs^-1 was used for the smallest, land-terminating glaciers (Ariebreen, Renardbreen, Scottbreen, and Werenskioldbreen). The choice of these values is discussed later. Sample radargrams can be found in figure 5 of Lapazaran et al. (2013), for two shallow profiles (max. thickness ∼35 and ∼75 m), in figure 2 of Jania et al. (2005) for a medium-thickness profile (∼200–290 m), and in figure 5 of Vasilenko et al. (2011) for a deep-bed profile (max. thickness ∼500 m).

As shown in Table 1, some glaciers were radio-echo sounded during several campaigns, involving time spans up to 7 years. To account for thickness changes during such periods, we applied the elevation change rate for Wedel Jarlsberg Land glaciers of −0.65 ± 0.08 m a−1 estimated by Nuth et al. (2010) differencing ICESat tracks as well as older topographic maps. We took as reference year for the correction that of the campaign with the largest amount of data collected for each particular glacier. For some glaciers we had to combine measurements taken in different seasons (summer versus spring), thus involving distinct snow-cover conditions. For instance, a crossover analysis of Ariebreen ice-thickness data revealed a 3-m difference between the data from the summer 2006 and the spring 2007 surveys. Based on snow thickness measurements made on the neighboring Hansbreen in April 2008 (Grabiec et al., 2011), this thickness difference was attributed to the spring snow cover, and the thickness data were adjusted accordingly.

Ice-thickness data were interpolated into a regular grid using ordinary kriging. Kriging is a geostatistical interpolation method in which the value at an unobserved location is predicted by weighted average of the values at surrounding locations, with weights selected according to a model fitted to a variogram that describes the spatial correlation of the data (Cressie, 1993). The use of kriging interpolation differs from a previous paper (Martín-Español et al., 2013), where we used the Gridfit routine available for MATLAB (D'Errico, 2006), choosing the triangular interpolation method and leaving the default values for the remaining parameters (smoothness = 1, regularizer = gradient, solver = normal). Gridfit is an approximation method (as opposed to interpolation) closely related to thin-plate spline approximation. When Gridfit is used to create soft surfaces, some of the field data can be ignored, resulting in large deviations of the interpolated surface from the individual measurements, largest in the roughest areas. This happens because Gridfit is a robust method to outliers and duplicates, which means that it ignores data that strongly deviate from the mean trend of the data set. However, this can constitute a source of error if such data are not detected and deleted a priori from the data set. Otherwise, the measurement error will be calculated from an improper data set, and the interpolation error estimate will also be inadequately influenced. Moreover, rough bedrock areas that have not been densely sampled will be seen as outliers and will be wrongly softened. In contrast to Gridfit, kriging is able to get reasonably smooth surfaces without losing reliability. In fact, it allows calibrating the maximum degree of deviation from the measurement data through a parameter called nugget (a selection of nugget different from zero means that kriging is no longer an exact interpolator). Kriging also provides a closer fit to the glacier physical conditions, allowing for anisotropy, through its ability to account for a higher autocorrelation in one direction than in another. A more detailed comparison between kriging and Gridfit, as concerns their suitability to calculate glacier volumes from GPR ice-thickness data, can be found in Martín-Español (2013, §3.6.1). It focuses on the capabilities of both methods to (1) work with sparse, error-containing and duplicate data; (2) deliver a smooth surface without losing reliability; (3) produce accurate ice-volume estimates on synthetic glaciers and ice caps; and (4) involve low interpolation errors.

We adjusted the 2007 glacier boundaries (Arendt et al., 2012; König et al., 2014) to the year of the GPR measurements using satellite imagery (mostly Aster and Landsat 7 images) to clip the gridded ice-thickness map. Glacier volume is the result of summing up the products, at each grid cell, of the cell thickness times the cell area. Ice volume estimations are affected by errors from various sources, which are discussed in the section Glacier Volume Error Estimates.

TABLE 1

Ground-penetrating radar (GPR) surveys, with listing of radar equipment used and total number of km of profiles done on each glacier. Ramac/GPR is manufactured by Malå Geoscience. VIRL2, VIRL6, and VIRL7 are GPR systems made in-house and described in Vasilenko et al. (2002), Berikashvili et al. (2006), and Vasilenko et al. (2011), respectively.

UNSURVEYED TRIBUTARY GLACIERS

We are interested in estimating the volumes of entire glacier basins. However, several of the investigated glaciers have tributaries that, due to time constraints or safety reasons, could not be surveyed. Although they represent a small percentage of the total volume, they have a severe effect on the volume-area ratio, and ignoring them would imply a bias in the volume estimate. If the thickness of the tributaries were interpolated from the closest GPR measurements on the main glacier trunk, plus zero-thickness data points at the walls of the tributary glacier, the volume of the tributaries would be systematically underestimated. To get a closer approximation to the volume of such tributary glaciers, we built an empirical function fitting the measured ice thickness along the center line of a set of 22 tributary glaciers in Svalbard for which GPR ice-thickness measurements were available. Such a function starts from zero thickness at the point of contact of the center line with the tributary glacier wall, and increases until the zone of confluence with the main glacier trunk. The measured ice thicknesses of each tributary, normalized by that of its deepest point (at the zone of confluence with the main trunk), as a function of the normalized distance to the tributary headwall, is shown in Figure 3, part a, together with its least-square fit to a third-degree polynomial.

This “tributary thickness function” provides an estimate of the expected ice thickness along the center line of the tributary glaciers, but we have to combine it with some estimate of the thickness variations along the transverse direction in order to calculate the volume of the tributary. Once the areas of different cross-sections along the tributary have been determined, the volume of the slice between two consecutive cross-sections is calculated as the product of their average area and the distance between them, and the total volume of the tributary is obtained as sum of the volumes of all of the slices.

There are different approaches in the literature to approximate the cross-section of valley glaciers as power law regressions (y = αx^β) and polynomial regressions (e.g., Doornkamp and King, 1971; Wheeler, 1984; James, 1996). However, there is not a definitive method to address this issue (Greenwood and Humphrey, 2002). We used a polynomial approach to estimate the area of the transverse cross-sections of the tributaries. We integrated, along the cross-section, the polynomial interpolating the ice thickness at the center of the cross-section (provided by the tributary thickness function) and the zero-thickness points at both extremes (the contact points with the valley walls). Having only three data points, there is a unique second-order polynomial passing through all of them. A fourth-degree polynomial has sometimes been found to give a better approximation to the cross-section of valley glaciers (Greenwood and Humphrey, 2002). However, there is a double-infinity set of fourth-order polynomials passing through three given data points. Accordingly, we tested a fourth-order polynomial excluding the nonsymmetric linear and cubic terms, so there is only one such polynomial passing through the three given points. The comparison of the calculated volumes of the echo-sounded tributaries (from the entire set of available ice-thicknesses for a given tributary plus zero-thickness boundary points at the contact of the tributary with its valley walls) with those estimated using the tributary thickness function and both second-order and fourth-order approaches for the cross-sectional areas showed a better agreement when the parabolic cross-section was used. The relative improvement in the volume estimate implied by using the parabolic approach versus using the fourth-order approach was, on average, 3% of the glacier volume, and increased to 5% for the smaller tributaries (volume < 0.3 km³). This is not surprising, since we are dealing with tributary glaciers, often small and more confined by the surrounding mountains as compared with the main trunk of the glacier. Their beds are therefore expected to be less flat than that of the main trunk, and hence the parabolic approach can be expected to work best. An exception to the parabolic cross-section was used in the present study: Recherchebreen has two large unsurveyed areas to the east, which are not typical small tributaries but, instead, rather large valley glaciers (see Fig. 4). Consequently, a fourth-degree polynomial approximation was used for calculating their cross-sections.

FIGURE 3.

(a) Normalized ice thickness H* as a function of the normalized distance x* to the headwall, for a collection of 22 Svalbard tributary glaciers. (b) Standard deviation σ_ΔH of the differences between the measured (GPR) and computed (tributary thickness function) ice thickness for each tributary glacier, and associated linear fit. (c) Tributary glacier cross section incremented/decremented by an amount equal to σ_ΔH.

The standard deviation (non-normalized) of the differences between the entire set of measured thicknesses at the tributaries and the corresponding ones estimated from the tributary thickness function is 32.4 m, which gives an idea of the error in thickness involved when using the tributary thickness function. We assess the error in the estimation of the volume of unsurveyed tributary glaciers as follows. For each tributary with available ice-thickness measurements, we compute the standard deviation σ_ΔH of the differences between the measured (GPR) and computed (tributary thickness function) non-normalized ice thickness along the center line of the tributary, and then plot such standard deviations versus the corresponding maximum ice thickness (Fig. 3, part b). A linear fit is then performed, which is used to estimate the error in ice thickness along the center line of the unsurveyed tributaries as a function of its maximum ice thickness, the latter being taken as the thickness measured at the closest GPR profile in the main glacier trunk. We still have to propagate this error to the transversal direction and to the volume computation. We do this as illustrated in Figure 3, part c: the tributary parabolic cross-section is incremented and decremented by an amount σ_ΔH, and the average of the associated volume changes along the entire tributary is taken as an estimate of the error in volume of the tributary.

RADIO-WAVE VELOCITY

When working with small glaciers and ice caps, the firn thickness distribution as well as the spatial distribution of cold and temperate ice are often unknown or only poorly resolved. Consequently, quite often the two-way travel time recorded by GPR is converted into ice thickness using a constant radio-wave velocity (RWV) for the entire glacier. However, it is well known that the RWV of glacier ice depends on many factors (most notably, air and water content fractions) and consequently changes substantially both in space and time. In any glacier, the firn layer has typically a RWV of 190–200 m µs, substantially larger than that of the ice (typically 160–170 m µs^-1). In polythermal glaciers such as those commonly found in Svalbard, the cold and temperate layers have different velocities, typically about 168–170 m µs^-1 for cold ice, depending on the air fraction, and 160–165 m µs^-1 for temperate ice, depending on the water content fraction. Temporal RWV variations are also common, mostly associated with the changing melting conditions and associated amount of water in the glacier body (Jania et al., 2005; Navarro and Eisen, 2010).

A proper glacier volume estimate would therefore require a good knowledge of the glacier internal structure (distribution and thickness of the firn, cold and temperate ice layers) and the ice and water fractions, which is seldom—if ever—available. With the aim of finding an optimal choice of constant RWV to be used for time-to-depth conversion, we carried out an experiment using information about the hydrothermal structure of Hansbreen, one of the best studied polythermal tidewater glaciers in Svalbard. Figure 5 displays the hydrothermal structure of Hansbreen along its center line, combining GPR data from Moore et al. (1999) and Pälli et al. (2003) with data from the authors. Moore's and Pälli's 50 MHz GPR data of May 1997 provided a better resolution for the firn layer thickness, while the 25 MHz data from the authors of April 2008 provided a more complete and updated picture of the cold layer thickness. The authors' 20 MHz GPR data collected from helicopter in April 2011 allowed determination of the ice thickness of the highly crevassed zone near the calving front. As precise data for the firn and cold-layer thickness are only available along the center line, their thickness was assumed to remain constant, for a given distance to the glacier front, across the glacier transversal section. Considering such a layered structure, our two-way travel time data collected during 2008–2011 were converted into ice thickness using RWVs of 190, 170, and 165 m µs for the firn, cold ice, and temperate ice layers, respectively. These choices were based on available layer velocities measured on Hansbreen by the authors using the common mid-point method (CMP); a sample CMP profile can be found in Figure 3 of Jania et al. (2005). The total resulting volume for Hansbreen was 10.74 ± 0.68 km³. We then performed volume estimates using constant RWVs for the entire glacier, between 165 and 170 m µs^-1, at 1 m µs^-1 intervals, and found that the constant RWV that produces the glacier volume closest to that computed using different velocities for distinct media is 166 m µs. This value is within the range of column-averaged RWVs measured on Hansbreen by the authors using the CMP method (Jania et al., 2005; Grabiec et al., 2012). The corresponding volume is 10.75 ± 0.73 km³, equivalent, within error bounds, to the volume given above. We thus selected this RWV as a glacierwide constant velocity for time-to-depth conversion of GPR data from the largest and thickest glaciers of our data set (Hansbreen, Paierlbreen, Recherchebreen, and Austre Torellbreen), all of which show a hydrothermal structure similar to that of Hansbreen and are of tidewater type except Recherchebreen. However, for the smaller and thinner glaciers of our data set (Ariebreen, Renardbreen, Scottbreen, and Werenskioldbreen), which show a smaller fraction of temperate ice (absent in the case of the smallest glacier, Ariebreen, which is entirely made of cold ice), it seems more realistic to use a higher RWV. We took 168 m µs^-1, on the basis of regional RWV measurements by the authors (Navarro et al., 2005; Jania et al., 2005).

FIGURE 4.

Ice-thickness map of the studied glaciers. The location of the two main panels shown is indicated in Figure 1. The green areas represent unsurveyed tributary glaciers whose volumes were estimated using the procedure described in the section Unsurveyed Tributary Glaciers. The inset illustrates, as an example, the ice-thickness distribution for a tributary glacier of Recherchebreen calculated using the mentioned procedure.

GLACIER VOLUME ERROR ESTIMATES

It is important that glacier volumes be accompanied by error estimates as accurate as possible. The error in glacier volume is computed, using standard error propagation techniques (e.g., Bevington and Robinson, 2002), from the errors of its main components. The detailed procedure is out of the scope of this paper, and we will restrict to present a conceptual approach, explaining how we take into account and combine together the main error sources. Further detail can be found in Martín-Español (2013, §3.7).

Given a digital ice-thickness model of a glacier, its volume can be computed as the sum of cell area times the average ice thickness for the cell. In a model, the grid cells have a predefined error-free area, and thus the error in volume would depend only on the errors in thickness at each individual grid point/cell. However, a digital ice-thickness model is constructed from a given boundary and a set of non-uniformly distributed ice-thickness data points obtained from GPR measurements along certain radar profiles. Furthermore, the real glacier boundary often cannot be accurately defined, because of, for example, the presence of snow patches covering the terrain surrounding the accumulation area, or the presence of debris cover on the ablation area, and this uncertainty in boundary delineation produces an error in volume that can be large. Once a given glacier boundary is assumed, it still has to be fitted using a certain grid (often regular), and this involves errors, both systematic and random, incurred by the mapping software used for creating the grid. Even if the systematic part (bias) is corrected, a random pixelation error still remains, which represents the error in volume at the boundary cells. On the other hand, determining the ice thickness at the grid points from the GPR measurement points (which are in turn affected by measurement errors) involves an interpolation error that can also be large. Estimating all of these errors, and their combined effect, is not straightforward. To help to do it, we will conceptually separate the errors in volume as those stemming from the uncertainty in boundary delineation, ϵ_VD, which we will estimate separately, and, once a given boundary is assumed, those incurred by all other computational aspects, ϵ_VC, including the pixelation error and those related to ice-thickness uncertainties, both measurement errors and interpolation errors.

FIGURE 5.

Polythermal structure along the center line of Hansbreen. In the ablation zone, a cold ice layer, with a maximum thickness of ca. 100 m some 5 km from the glacier front, overlies a temperate ice layer. The accumulation zone has a firn layer with average thickness of 15 m. The ice underlying the firn layer is temperate.

Regarding the uncertainty in boundary delineation, glacier outlines are often defined from aerial photos or satellite images taken at the end of the summer, when the snow cover on the terrain surrounding the glacier is expected to be minimal. But, even so, on the valley walls surrounding the accumulation area often remain snow patches that make difficult the proper definition of the glacier boundary. Similarly, on the ablation area, near the contact with the glacier valley walls and also near the snout of land-terminating glaciers, the presence of debris makes difficult a proper interpretation of the glacier limits. These errors are strongly glacier-dependent (e.g., they depend on glacier orientation and altitude of surrounding mountains, and associated shadowing, and on the slope and petrology of glacier valley walls, among other factors) and their estimation is subject to strong uncertainties, even in the case of boundaries delineated from ground-based measurements such as GPS positioning of the glacier outlines. Consequently, we have taken as uncertainty in boundary delineation the 8% error in area estimated for the most recent Svalbard glacier inventory (König et al., 2014). This error in area implies a corresponding error in volume, which could at a first glance be thought to be small, because the uncertainty in boundary delineation occurs where the ice thickness is expected to be smallest, namely near the zero-thickness points at the contact with the glacier valley walls and at the glacier front of land-terminating glaciers. However, the relatively scarce coverage of the glacier surface by the radar profiles, especially near the glacier walls, implies that an improper definition of the glacier limits has an impact on the ice-thickness interpolation that extends well into the glacier. It is shown in Martín-Español (2013, §3.7.3) that the error in volume incurred by the uncertainty in boundary delineation can be roughly approximated, for a glacier with a radar profile net such as those of the glaciers included in this study, by the product of the error in area (in our case, the 8% error mentioned above) and one-half of the average thickness of the glacier.

Once a given boundary delineation is assumed, and supposing that it is known how the planar area for a given two-dimensional boundary is computed by the software used to construct the grid, it is possible to correct for the bias in area introduced by the mapping software and the magnitude of the random pixelation error, ϵVPix, that remains after the removal of the systematic error. Both of them become smaller with the reduction of the grid size. The pixel size of the satellite image, or the horizontal positioning error, in the case of ground-based boundary delineations, play the role of a lower bound in the computation of the pixelation error, so the grid size should be chosen accordingly to avoid unnecessary computational load (the grid size should not be lower than the pixel size or the error in horizontal positioning). After correcting for biases, the magnitude of the pixelation error is usually negligible as compared with the uncertainty in boundary delineation discussed above, as is the case with our glacier collection for Wedel Jarlsberg Land.

Concerning the errors in thickness, a first step is to perform a crossover analysis aimed to check for possible discrepancies between the ice-thickness data at intersecting points between profiles. Such discrepancies can be caused by inaccurate positioning of the data, errors in the picking of the bed reflections, or lack of (or improper) migration of any of the intervening profiles. The latter is known to have a strong impact on the thickness estimate for highly sloping beds. The apparent change in thickness due to migration can be approximated by the cosine of the angle between the surface and bed slopes in the direction of profiling. This angle often attains its highest values for profiles perpendicular and close to the valley walls. At high angles (e.g., 30–45°), the ice thickness for the migrated profiles can increase up to 25–40% with respect to that of non-migrated profiles. Although these are extreme cases, our transverse profiles for Ariebreen experienced changes up to 30% at locations close to its steep valley walls. In general, however, the thickness changes implied by migration were of the order of 3–4%. Though the bed slopes close to the valley walls of the studied glaciers were often large, the GPR profiles perpendicular to them had to be interrupted in many cases at some distance from the side wall because of either large surface slopes (often involving deep fresh snow) or presence of lateral crevasses.

Assuming properly migrated data, if the intersecting profiles correspond to a single radar campaign, the discrepancies at the crossovers should be of the order of the errors in ice thickness discussed below. If, however, they correspond to different radar campaigns (as happened in many of our studied glaciers), these discrepancies will show a systematic pattern (bias) that should be corrected as discussed earlier. In what follows, we assume that such biases have been corrected for. Even if corrected for such biases, our migrated GPR data still showed high crossover differences for crossovers involving radar profiles parallel and close to steep valley side-walls. This setting implies that the first arrivals from the bed can correspond to off-nadir reflections from the valley walls. 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 it often fails for profiles parallel to them, and hence a large difference appears at crossovers close to the valley walls. This could be corrected by applying three-dimensional migration (e.g., Moran et al., 2000). However, our net of radar profiles was not dense enough to allow for three-dimensional migration. Consequently, to avoid these problems we discarded the radar profiles parallel and close to the side-walls. This resulted in overall crossover differences of the order of the errors in ice thickness discussed below.

The errors in ice thickness arise from different sources. First, we have measurement errors, which embrace GPR errors and GPS errors. GPR measurements are done at particular points, and their output are two-way travel times, which are converted to ice thickness multiplying them by half the RWV. Consequently, there is a point-dependent GPR error that depends on both the error in the assumed RWV and the errors in timing. The latter, in turn, include the errors associated with the time discretization of the radar signal (though these are actually superseded by the range or vertical resolution of the radar) and the picking errors incurred when improperly picking the bed reflections from the radargram. The latter are interpretation errors that are minimized when working with properly migrated radargrams, as has been done with our GPR data set. Following Navarro and Eisen (2010), we estimated the GPR error assuming that the error in two-way travel time is given by the inverse of the central frequency of the radar system and that the error in RWV is typically of the order of 2%. But, additionally, the location of the measurement points is affected by a horizontal positioning error that induces, for sloping surface/bed, an error in thickness (GPS error). We used, for radar-trace positioning, stand-alone GPS systems, for which we can assume a typical horizontal accuracy of 5 m (Wing et al., 2005). To estimate the thickness change that could be expected from such an error in horizontal positioning we need to know how much ice-thickness variability exists in 5 m around the measurement points. This variability could theoretically be measured as the standard deviation of the differences in ice thickness between pairs of data points 5 m apart from each other. As it is not common to have available such measurements, we estimate the variability from the interpolated ice-thickness map. Following Martín-Español (2013, §3.7.2), we evaluate the standard deviation of the difference in ice thickness between the measurement point and at least eight surrounding points situated 5 m away on the digital ice-thickness map, and take this value as the error in thickness introduced by the uncertainty in horizontal GPS positioning. This error increases with the difference in slope between glacier surface and bedrock.

Second, the ice thickness at the grid points are calculated, from the thickness data at the GPR measurement points, by means of an interpolation method that involves an interpolation error. The latter is an error of method, which will produce nonzero errors for the interpolated ice thickness at the grid cells even in the hypothetical case of error measurement-free GPR-retrieved ice-thickness data. A proper estimate of the interpolation error should relate well to the quantity and the distribution of field data available. Our study implements a novel method to account for interpolation error that has been designed to work with data distributions typical of GPR surveys, which show a high density of data along particular directions (profiles), often intersecting, but lack data in the areas between profiles. Our method aims to calculate 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 neighbor. We first construct a function relating the interpolation error at any given point with the distance to the nearest GPR-measured data point (distance-error function, DEF). Then, we calculate a glacier-wide interpolation error by averaging the interpolation errors computed by the DEF at each grid point. This method has the added benefit that it also provides a way to detect any possible bias involved in the interpolation process. Further details about this technique can be found in Martín-Español (2013, §3.7.2).

Each GPR measurement point has an associated measurement error that has been discussed earlier. This error will propagate to the thickness estimates at the grid cells through the interpolation algorithm used. Consequently, we estimated the propagation of the measurement error using the same neighbor distribution and weights as we did for the interpolation of ice thickness.

The combined error in ice thickness at a given grid cell is then calculated as root of the sum of squares of the (propagated) measurement error and the interpolation error at the cell.

As far as how to combine all of the above errors to estimate the error in glacier volume, we first note that the boundary delineation error ϵ_VD and the computation error ϵ_VC are nearly independent, so they can be combined as root of the sum of squares. The part of the computation error associated to the pixelation error is negligible (as compared to the other error sources) in our study cases. The ice-thickness-related error is thus the dominant contributor to the computation error ϵ_VC. The ice-thickness estimates at the grid cells are not independent of each other. Martín-Español (2013, §3.7.3) showed that the error in volume stemming from errors in thickness at the grid cells can be approximated by the product of the glacier area and a glacier-wide ice-thickness error (computed as standard deviation of the random variable error in thickness at the grid cells), divided by the square root of the maximum number of grid nodes within the glacier that are separated by a distance larger than the range of the kriging variogram, where we can assume spatial independence.