We present ground-penetrating radar (GPR)—based volume calculations, with associated error estimates, for eight glaciers on Wedel Jarlsberg Land, southwestern Spitsbergen, Svalbard, and compare them with those obtained from volume-area scaling relationships. The volume estimates are based upon GPR ice-thickness data collected during the period 2004–2013. The total area and volume of the ensemble are 502.91 ± 18.60 km^{2} and 91.91 ± 3.12 km^{3}, respectively. The individual areas, volumes, and average ice thickness lie within 0.37–140.99 km^{2}, 0.01–31.98 km^{3}, and 28–227 m, respectively, with a maximum recorded ice thickness of 619 ± 13 m on Austre Torellbreen. To estimate the ice volume of unsurveyed tributary glaciers, we combine polynomial cross-sections with a function providing the best fit to the measured ice thickness along the center line of a collection of 22 surveyed tributaries. For the time-to-depth conversion of GPR data, we test the use of a glacierwide constant radio-wave velocity chosen on the basis of local or regional common midpoint measurements, versus the use of distinct velocities for the firn, cold ice, and temperate ice layers, concluding that the corresponding volume calculations agree with each other within their error bounds.

## Introduction

The downward revision by Shepherd et al. (2012) of some previous estimates of the current contribution of the ice sheets to sea-level rise has renewed the interest on the contribution by glaciers and ice caps (GICs)—though, as pointed out by Hanna et al. (2013), the most recent estimates for the latter (e.g., Jacob et al., 2012; Gardner et al., 2013) also show values lower than previous estimates (e.g., Cogley, 2012). The potential contribution of GICs to sea-level rise depends on their total volume, and therefore the estimation of the total ice volume stored by GICs has recently received much attention (Radić and Hock, 2010; Huss and Farinotti, 2012; Grinsted, 2013). Because of the scarcity of ice-thickness measurements by methods such as ground-penetrating radar (GPR), seismic soundings, or deep ice drilling, the global estimates of the volume of GICs are based either on volume-area (V-A) scaling relationships (Radić and Hock, 2010; Grinsted, 2013) or on physically based simple models relating ice-thickness distribution to glacier geometry, mass balance, and dynamics (Huss and Farinotti, 2012). These two global methods still rely on the availability of accurate volume estimates from ice-thickness measurements, which are used as input data for the V-A relationships, or for tuning the physically based relationships. Volume-area and volume-length scaling relationships are also a usual tool for regional and global projections of glacier wastage, and corresponding sea-level rise, in response to climate scenarios (Radić et al., 2013).

The main aim of the present paper is to contribute to the enrichment of the scarce global data base of available glacier volumes retrieved from observations of ice thickness. We here present volume calculations, with associated error estimates, for eight glaciers on Wedel Jarlsberg Land, southwestern Spitsbergen, Svalbard (Fig. 1). At a regional scale, these new volume estimates add to the available set of GPR-based volume estimates of individual glaciers on Svalbard (Dowdeswell et al., 2008; Pettersson et al., 2011; Lapazaran et al., 2013; Saintenoy et al., 2013; Martín-Español et al., 2013). These estimates have allowed a new V-A relationship specific for Svalbard glaciers (Martín-Español, 2013) to be derived that is more reliable than those currently available (Macheret and Zhuravlev, 1982; Hagen et al., 1993), which were based on the limited and low-accuracy ice-thickness estimates from airborne radio-echo soundings in the late 1970s and early 1980s (Macheret, 1981; Macheret and Zhuravlev, 1982; Dowdeswell et al., 1984a, 1984b). The interest of V-A scaling relations characterizing individual glacier shapes, slopes, and sizes, as well as plausible glaciological conditions (such as steady-state, or sustained retreat or advance), has recently been pointed out by Adhikari and Marshall (2012), and thus regional V-A relations have become a matter of particular interest.

In addition to the interest of the glacier volumes presented here, there are many studies which require as input data the spatial distribution of ice thickness, as provided by the GPR data reported in this paper. Most relevant are the models, including a dynamical component, aimed to project the future contribution of glaciers and ice caps to sea-level rise (e.g., Radić et al., 2007, 2008) or the changes in the hydrology of glacierized mountain regions (e.g., Huss et al., 2008). A second aim of this paper is therefore to make available such ice-thickness distribution data for some Svalbard glaciers.

## Field Data

The GPR ice-thickness data used in this study correspond to eight glaciers located on Wedel Jarlsberg Land, southwestern Spitsbergen. Their locations, and their outlines in 2007 according to the Randolph Glacier Inventory (Arendt et al., 2012), are shown in Figure 1. The data were collected during several field campaigns between 2004 and 2013, using different radar equipment and frequencies, which are summarized in Table 1. The majority of the field campaigns were carried out in spring, before the onset of melting (about 95% of the profiles). In general, basal reflections were easily observed in all measurements done with VIRL6 and VIRL7 20 MHz radar systems, while those done using 25 MHz Ramac/GPR often did not allow clear bed detection for ice thickness above ∼300 m. These glaciers have very different characteristics in terms of size, morphology, and dynamics. Their surface areas span from 0.37 km^{2} up to 141 km^{2}. Five of them are land-terminating, while the remaining three are tidewater. They all have thinned since 1990 (Nuth et al., 2010). One of them, Paierlbreen, experienced a surge during 1992–1996 (Błaszczyk et al., 2009).

Figure 2 shows the distribution of the GPR profiles over the glacier surfaces. Although most glaciers are covered by a sufficiently dense net of profiles, there are some areas that were not surveyed due to time limitations or the presence of crevasses. Because our goal is to calculate the ice volume of entire glacier basins, we circumvent this problem by providing approximate ice-thickness data for these unsurveyed areas, as described in the section Unsurveyed Tributary Glaciers.

## Methods

### 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^{3}). 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.

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^{3}. 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^{3}, 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).

### 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, ϵ

*, including the pixelation error and those related to ice-thickness uncertainties, both measurement errors and interpolation errors.*

_{VC}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.

## Results and Discussion

The computed areas and volumes for the studied glaciers, together with their error estimates, are shown in Table 2, which also includes the maximum and average thickness for each glacier. The associated ice-thickness maps are displayed in Figure 3. The total area and volume of the ensemble are 502.91 ± 18.60 km^{2} and 91.91 ± 3.12 km^{3}, respectively. The individual areas, volumes, and average thickness show large spans, of 0.37–140.99 km^{2}, 0.01–31.98 km^{3}, and 28–227 m, respectively. The maximum ice thickness, reaching 619 ± 13 m, is found in the northernmost part of Austre Torellbreen, near its ice divide with Vestre Torellbreen, within the Amundsenisen ice field.

The volume of Hansbreen presented here is about 1 km^{3} (9%) larger than that recently published by Grabiec et al. (2012), who used the same set of radar profiles. The main reason is that the outlines of this glacier in the Randolph Glacier Inventory (used here) and in Grabiec et al. (2012) differ substantially at the ice divide between Paierlbreen and Hansbreen, implying a difference in the surface area of Hansbreen, which is some 8 km^{2} larger in the Randolph Glacier Inventory. Additionally, Grabiec et al. (2012) used a lower RWV (164 m µs^{-1}) for the time-to-depth conversion, producing smaller ice thicknesses, which also contribute to make the glacier volume smaller.

The relative errors in volume for the individual glaciers are all within 4–8%. The two main components of the error in volume (uncertainty in boundary delineation, ϵ_{VD}, and computational errors, ϵ_{VC}) contribute similar shares to the total error in volume. In the case of Recherchebreen, the large share of the computational error is due to the large size of the unsurveyed areas, for which the tributary thickness function and polynomial cross-sections were used to approximate their volume.

As an alternative to using the tributary thickness function, combined with second- or fourth-degree polynomial cross-sections, to estimate the volume of the unsurveyed tributary glaciers, we tested the use of a general V-A scaling relationship derived by Grinsted (2013) using the Randolph Glacier Inventory. We applied both methods to our set of 22 surveyed tributary glaciers and compared their results with the volumes calculated from the GPR ice-thickness data. Using the V-A scaling resulted in volumes lower by 50% on average than the GPR-derived volumes (i.e., a significant underestimation), with standard deviation of 66%, as compared to volumes lower by only 2% on average (indicating no significant bias), with standard deviation of 24%, when using the tributary function approach. The likely reason for the poor behavior of V-A scaling for estimating the volumes of tributary glaciers is that the scaling relationships are derived from volumes of entire glacier basins, and consequently they do not fit well to the case of tributary glaciers, which show particular aspect and thickness ratios. The use of the tributary thickness function, combined with second-degree (for small tributaries) or fourth-degree (for large tributaries) polynomial cross-sections, is therefore a clear improvement for the estimation of the volumes of the unsurveyed tributary glaciers.

The similar volumes obtained for Hansbreen when using, for the time-to-depth conversion, a proper choice of constant “average” RWV for the entire glacier or different RWVs for the distinct media (firn, cold ice, temperate ice) illustrate that, when using the simpler approach of constant RWV, the ice thicknesses overestimated in the firn and cold ice layers (where the real RWVs are higher than the average used as an approximation) are nearly compensated by the underestimation occurring at the temperate ice layer (where the real RWV is slightly lower than the average). A given relative error in RWV produces relative errors in ice thickness and volume of the same magnitude. Consequently, the volume estimate is not extremely sensitive (compared to other error sources) to the choice of a constant RWV for the entire glacier, provided that the latter is based on local or regional RWV values measured in similar seasons to those of the ice-thickness measurements, to minimize temporal variations of the RWV.

## TABLE 2

Calculated area and volume for the studied glaciers, with their sums and error estimates. The contributions to the error in volume associated with uncertainty in boundary delineation (ϵ_{VD}) and with computational aspects (ϵ_{VC}, dominated by thickness-related errors) are detailed, as well as the percentage of each glacier's volume corresponding to tributary glaciers. The average and maximum thickness of the glaciers are also given, together with their error estimates.

Different volume-area scaling relationships exist in the literature. We here analyze how these relationships perform in this particular region, with the awareness that V-A relations are expected to produce proper volume estimates only when applied to large sets of glaciers, and not to individual glaciers or small sets of glaciers. Figure 6 shows the relative error (as percentage) of the volume estimates produced by several scaling relationships when applied to each of our measured glaciers. In contrast with the volume results from V-A scaling relationships obtained by Martín-Español et al. (2013) for glaciers in western Nordenskiöld Land, central Spitsbergen, we observe, for Wedel Jarlsberg Land glaciers, an almost systematic underestimation of the glacier volumes, except in the case of the Radić and Hock (2010) relationship, which overestimates the total volume of the ensemble by 12%. This particular scaling relationship has been suggested to overestimate the volumes (Grinsted, 2013). The underestimate of Wedel Jarlsberg Land glacier volumes by most of the V-A relationships is likely due to the fact that these glaciers are generally thicker than those in western Nordenskiöld Land, and the share of tidewater glaciers in the Wedel Jarlsberg Land sample is far larger than in the western Nordenskiöld Land study, so the real volume/area ratio is expected to be higher. The best overall result is provided by the relationship derived by Hagen et al. (1993), which underestimates by 4% the measured volume of the ensemble, whereas the worst performance is shown by the relationship of Macheret and Zhuravlev (1982), which underestimates the volume of the ensemble by 34%. Dowdeswell et al. (1984a) pointed out that the results from airborne Soviet radio-echo soundings in the early 1980s using 440 and 620 MHz radar systems (Macheret, 1981; Macheret and Zhuravlev, 1982) had underestimated the ice thickness by 33–50%, and consequently their empirical relationship based on the measured volumes is biased in a similar way. The other scaling approaches that underestimate the volume of the ensemble are those of Chen and Ohmura (1990), by 19%; Bahr et al. (1997), by 21%; and Grinsted (2013), by 26% when using the relationship derived from the GLIMS Inventory, and by 10% when using the relationship derived from the Randolph Glacier Inventory. Concerning particular glaciers, the volume of Paierlbreen is overestimated by all V-A relations, by 16% on average (the largest average overestimate in our set of glaciers). We attribute this to the recent surge, in the 1990s, of Paierlbreen, because a surging glacier, in post-surge phase, has a lower volume/area ratio due to the frontal advance and thinning after the surge event (Cuffey and Paterson, 2010, chapter 12).

The glacier volume data presented here have been used, together with additional data from the authors or from other authors for different Svalbard regions, to build a new regional V-A scaling relationship for Svalbard glaciers that has been used to calculate the volume of the entire population of Svalbard glaciers, estimated to be 7504 ± 312 km^{3}, or 19 ± 1 mm sea-level equivalent (Martín-Español, 2013).

To our knowledge, the hydrothermal structure of the glaciers studied in this paper has only been discussed in the literature for Hansbreen (e.g., Dowdeswell et al., 1984a, 1984b; Macheret et al., 1993; Jania, 1996; Jania et al., 1996; Moore et al., 1999; Pälli et al., 2003; Grabiec et al., 2012), Recherchebreen (e.g., Macheret and Zhuravlev, 1982; Dowdeswell et al., 1984a, 1984b), Werenskioldbreen (e.g., Macheret and Zhuravlev, 1982; Dowdeswell et al., 1984a, 1984b; Pälli et al., 2003), and Ariebreen (Lapazaran et al., 2013), all of them polythermal except the latter, which is entirely made of cold ice. On the basis of the presence or absence of scattering in the radargrams as an indicator of temperate or cold glacier ice (Navarro and Eisen, 2010), we here interpret that all other glaciers in our data set (Austre Torellbreen, Paierlbreen, Renardbreen, and Scottbreen) have a polythermal structure.

## Conclusions

We summarize as follows the main conclusions of our analysis:

(1)The total area and volume of the studied glaciers are 502.91 ± 18.60 km

^{2}and 91.91 ± 3.12 km^{3}, respectively. The individual areas, volumes, and average thickness lie within 0.37–140.99 km^{2}, 0.01–31.98 km^{3}, and 28–227 m, respectively. The maximum measured ice thickness (619 ± 13 m) was recorded on Austre Torellbreen.(2)The error estimates for the studied glaciers represent 4–8% of their volume. The errors in volume due to the uncertainty in boundary delineation are of the same order of magnitude as the computational errors in volume associated with the measurement and interpolation errors.

(3)The use of a tributary thickness function describing a typical longitudinal ice-thickness profile of the tributary glaciers in this region, combined with polynomial cross-sections, produces a much better approach to the estimation of the volume of the unsurveyed tributary glaciers than the use of a general volume-area relationship.

(4)The use for the time-to-depth conversion of GPR data of a glacier-wide constant radio-wave velocity based on local or regional common midpoint measurements, versus the use of different velocities for the different materials (firn, temperate ice, cold ice), produces glacier volume estimates that agree with each other within their error bounds.

(5)All V-A relationships underestimate the total volume of our ensemble of glaciers, except the Radić and Hock (2010) relationship, which overestimates it by 12%. The best result is shown by the Hagen et al. (1993) relationship (4% underestimate), whereas the worst is shown by the Macheret and Zhuravlev (1982) relationship (34% underestimate).

## Acknowledgments

This research was supported by grant EUI2009-04096 (PolarCLIMATE-SvalGlac) from the Spanish EuroResearch Programme; grants CGL2005-05483, CTm^{2}008-05878, and CTm^{2}011-28980 from the Spanish National Plan for R&D; grant NCBiR/PolarCLIMATE-2009/2-2/2010 from the Polish National Centre for R&D; grants IPY/269/2006 and N N306 094939 from the Polish Ministry of Science and Higher Education; Polish-Norwegian funding through the AWAKE (PNRF-22-AI-1/07) project; and the ice2sea programme from the European Union 7th Framework Programme, grant number 226375, ice2sea contribution number 47. Satellite imagery used in this paper were available from ASTER © METI and NASA (2005), for August 2004 and July 2006, and from SPOT-5 © CNES (2008), for August 2008, distribution Spot Image S.A., all rights reserved. The paper was greatly improved by the comments and suggestions by Reinhard Drews and two anonymous reviewers. We also thank the associate editor, Jacob Clement Yde, for his efforts in improving the quality of the paper.

## References Cited

*Geophysical Research Letters*, 39: L16505, doi: http://dx.doi.org/10.1029/2012GL052712. Google Scholar

*Journal of Geophysical Research*, 102(B9): 20355–20362. Google Scholar

*Radiotekhnika*, 9: 52–57 (in Russian). Google Scholar

*Data Reduction and Error Analysis for the Physical Sciences*. Third edition. New York: McGraw-Hill, 320 pp. Google Scholar

*Polish Polar Research*, 30(2): 85–142. Google Scholar

*In*H. Lang , and A. Musy (eds.),

*Hydrology in Mountainous Regions (Proceeding of Two Lausanne Symposia, August, 1990).*IAHS Publication No. 193: 127–135. Google Scholar

*In*A. Henderson-Sellers , and K. McGuffie (eds.),

*The Future of the World's Climate*. Second edition. Amsterdam: Elsevier, 197–222. Google Scholar

*Statistics for Spatial Data*. Revised edition. Wiley Series in Probability and Mathematical Statistics. New York: Wiley. ISBN 0471002550. Google Scholar

*The Physics of Glaciers*. Fourth edition. Burlington: Elsevier, 704 pp. Google Scholar

*Numerical Analysis in Geomorphology: an Introduction*. London: Edward Arnold, 372 pp. Google Scholar

*Journal of Glaciology*, 30(104): 295–314. Google Scholar

*Norsk Polarinstitutt Skrifter,*182: 41 pp. Google Scholar

*Journal of Geophysical Research*, 113: F03022, doi: http://dx.doi.org/10.1029/2007JF000905. Google Scholar

*Science*, 340: 852–857. Google Scholar

*Polish Polar Research*, 32(4): 393–421. Google Scholar

*Polish Polar Research*, 33(2): 112–138. Google Scholar

*Computing Science and Statistics*, 34: 452–460. Google Scholar

*The Cryosphere*, 7(1): 141–151. Google Scholar

*Norsk Polarinstitutt Meddelelser,*129: 160 pp. Google Scholar

*Nature*, 498: 51–59. Google Scholar

*Journal of Geophysical Research*, 117(F4): doi: http://dx.doi.org/10.1029/2012JF002523. Google Scholar

*Hydrological Processes*, 22(19): 3888–3902. Google Scholar

*Nature*, 482: 514–518. Google Scholar

*Earth Surface Processes and Landforms*, 21(5): 413–432. Google Scholar

*Polar Research*, 15(1): 53–66. Google Scholar

*Annals of Glaciology*, 42(1): 125–134. Google Scholar

*In*J. S. Kargel , G. J. Leonard , M. P. Bishop , A. Kääb , and B. H. Raup (eds.),

*Global Land Ice Measurements from Space*. London: Springer-Praxis, 400 pp. Google Scholar

*Razvedka i Okhrana Nedr,*21: 6–11 (in Russian). Google Scholar

*Polar Research*, 32: 11068, doi: http://dx.doi.org/10.3402/polar.v32i0.11068. Google Scholar

*Annals of Glaciology*, 2: 45–51. Google Scholar

*Journal of Glaciology*, 28(99): 295–314. Google Scholar

*Journal of Glaciology*, 39(132): 373–384. Google Scholar

*Annals of Glaciology,*54(64): 211–217, doi: http://dx.doi.org/10.3189/2013AoG64A109. Google Scholar

*Journal of Glaciology*, 45(151): 524–532. Google Scholar

*Journal of Glaciology*, 46(153): 274–286. Google Scholar

*In*P. Pellikka , and W. G. Rees (eds.),

*Remote Sensing of Glaciers—Techniques for Topographic, Spatial and Thematic Mapping*. Leiden: CRC Press, 195–229. Google Scholar

*Annals of Glaciology*, 42(1): 158–162. Google Scholar

*Journal of Geophysical Research*, 115(F1): doi: http://dx.doi.org/10.1029/2008JF001223. Google Scholar

*Polar Research*, 22(2): 355–371. Google Scholar

*Geografiska Annaler*, Series A, Physical Geography, 93(4): 311–322. Google Scholar

*Journal of Geophysical Research*, 115(F1): doi: http://dx.doi.org/10.1029/2009JF001373. Google Scholar

*Annals of Glaciology*, 46: 234–240. Google Scholar

*Journal of Glaciology*, 54(187): 601–612. Google Scholar

*Climate Dynamics*, 42: doi: http://dx.doi.org/10.1007/s00382-013-1719-7. Google Scholar

*Near Surface Geophysics*, 11: doi: http://dx.doi.org/10.3997/1873-0604.2012040. Google Scholar

*Science*, 338: 1183–1189. Google Scholar

*Geophysics*, 44(6): 1041–1063. Google Scholar

*Bulletin of the Royal Society of New Zealand*, 35: 611–618. Google Scholar

*Journal of Glaciology*, 57(206): 1113–1118. Google Scholar

*Earth Surface Processes and Landforms*, 9(4): 391–394. Google Scholar

*Journal of Forestry*, 103(4): 169–173. Google Scholar

*Seismic Data Analysis: Processing, Inversion and Interpretation of Seismic Data.*Volumes 1 and 2. Tulsa, Oklahoma: Society of Exploration Geophysicists, 2027 pp. Google Scholar