## Introduction

Mass loss from glaciers is currently one of the largest contributors to sea-level rise, with a 27% share of the sum of the total contributions over the period 1993–2010 (Reference StockerStocker and others, 2013). Moreover, modelling studies show that the contribution to sea-level rise by glaciers will continue to be important during the 21st century. A recent multi-model study suggests a sea-level rise of 155 ± 141 mm and 216 ± 144 mm over the period 2006–2100, for climate scenarios RCP4.5 and RCP8.5, respectively, involving a reduction of the current global glacier volume by 29% and 41% (Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others, 2014). How much and for how long glaciers will continue to be an important contributor to sea-level rise depends on their total volume, the estimate of which has therefore recently received much attention (Reference Radić and HockRadić and Hock, 2010; Reference Huss and FarinottiHuss and Farinotti, 2012; Reference Marzeion, Jarosch and HoferMarzeion and others, 2012; Reference GrinstedGrinsted, 2013; Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others, 2014).

Due to logistic difficulties, glacier volumes determined from ice thickness measurements by methods such as ground-penetrating radar (GPR), seismic soundings or deep ice drilling are available for an extremely small fraction of the total population of glaciers. Adding our 25 glacier volumes reported here to the 337 glacier volumes available in the updated version of the catalogue by Reference Cogley, Henderson-Sellers and McGuffieCogley (2012) gives a total of 362 glacier volumes out of a global population of just under 200 000 glaciers in the recent Randolph Glacier Inventory (RGI) V3.2 (Reference PfefferPfeffer and others, 2014), i.e. <0.2%. This small representation highlights the importance of exploring indirect methods to approximate glacier volumes from other known parameters. There are two main indirect methods of inverting for volume from surface properties. Volume–area (*V*–*A*) scaling relationships follow from a dimensional analysis of the driving equations of glacier dynamics (Reference Bahr, Meier and PeckhamBahr and others, 1997), while physically based numerical inversion methods are based on numerical inversions relating ice thickness distributions to glacier geometry and dynamics, mass balance and thinning rates (e.g. Reference McNabbMcNabb and others, 2012). Both methods suffer from the same ill-posed nature of the inversion problem, as all boundary conditions are specified at the glacier surface (Reference Bahr, Pfeffer and KaserBahr and others, 2014).

Scaling relationships are easily implemented and have proved useful when dealing with large ensembles of glaciers, given their ability to characterize global and regional ice volumes (Reference Radić and HockRadić and Hock, 2010; Reference Marzeion, Jarosch and HoferMarzeion and others, 2012; Reference GrinstedGrinsted, 2013; Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others, 2014). The information required for this method includes glacier area and sometimes other glaciological parameters such as length or characteristic shape (width to length ratio), which can be extensively determined from satellite images and digital elevation models. Currently this information is freely available in global inventories such as the RGI, or can be easily derived from them. In contrast to scaling methods, physically based methods allow the calculation of ice-thickness distribution and subglacial topography. Because of the large amount of information usually required (e.g. surface topography, thinning rates, surface mass balance and surface velocity field (Reference McNabbMcNabb and others, 2012)), physically based methods are often restricted to the analysis of unique or small sets of glaciers (e.g. Reference Farinotti, Huss, Bauder and FunkFarinotti and others, 2009). However, recent applications have dealt with much larger sets of glaciers, at both regional (Reference FreyFrey and others, 2013; Reference Huss and FarinottiHuss and Farinotti, 2014) and global (Reference Huss and FarinottiHuss and Farinotti, 2012) scales, though the reduced set of input data required comes at the expense of general assumptions about variables such as surface mass balance and calving flux (Reference Huss and FarinottiHuss and Farinotti, 2012), whose parameterizations have been calibrated using small empirical datasets.

Svalbard is a highly glacierized archipelago situated in the Atlantic sector of the Arctic (76–81° N, 10–33° E; Fig. 1), which is a region highly vulnerable to climate change (ACIA, 2005). A recent study by Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others (2014), based on multiple climate models, projected mass losses from Svalbard glaciers, over the period 2006–2100, of 12.41 and 15.81 mm sea-level equivalent (SLE) for climate scenarios RCP4.5 and RCP8.5, respectively, representing reductions of 55% and 70% of the current total volume of Svalbard glaciers.

The ice masses of Svalbard cover ∼33 922 km^{2} (Reference PfefferPfeffer and others, 2014), putting this among the largest glacierized areas in the Arctic. The volume of the entire population of Svalbard glaciers has recently been derived, as part of worldwide studies, from global scaling relationships (Reference Marzeion, Jarosch and HoferMarzeion and others, 2012; Reference GrinstedGrinsted, 2013; Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others, 2014) and physically based methods as in Reference Huss and FarinottiHuss and Farinotti (2012), revealing high variance across different estimates (Table 1). Older estimates derived from Svalbard-specific scaling relationships (Reference YuYa and ZhuravlevMacheret and Zhuravlev, 1982; Reference Hagen, Liestøl, Roland and JørgensenHagen and others, 1993) led to more consistent results. More recent scaling relationships are mostly based on echo sounding flights carried out in the 1970s and 1980s using old radar and positioning systems and mostly consisting of a single profile along the centre line of each glacier (Reference YuYa and ZhuravlevMacheret and Zhuravlev, 1982; Reference Dowdeswell, Drewry, Liestøl and OrheimDowdeswell and others, 1984). The corresponding volume estimates, performed assuming parabolic cross-sections, are expected to have substantial errors. Currently, there is a greater amount of radio-echo sounding data that covers the majority of the surface of the surveyed glaciers. This, combined with improved positioning and recording systems, leads to more accurate volume estimates.

In this study we bring together all the available GPR datasets for Svalbard and calculate their associated glacier volumes, which we then use to calibrate new scaling relationships specific to Svalbard glaciers. We then apply these relationships to the most up-to-date glacier inventory, the RGI V3.2 (Reference PfefferPfeffer and others, 2014), to estimate the total volume of Svalbard glaciers and subsequently their total *potential* contribution to sea-level rise. Austfonna and Vestfonna are two large ice caps on Nordaustlandet (Fig. 1), which contain much of the glacier volume of Svalbard and have been extensively radio-echo sounded, thus allowing an accurate volume estimate. Since, on the other hand, the number of echo sounded ice caps on Svalbard is insufficient to build a separate *V*–*A* relationship for ice caps, Austfonna and Vestfonna are treated separately, simply adding their GPR-based volume estimates to the *V*–*A* relationship-derived volume estimate for the rest of Svalbard. Jan Mayen is a small (∼377 km^{2}) island located some distance (70°56′ N, 8°32′ W) from the Svalbard archipelago but is grouped together with it in the RGI. Consequently, we have included its 48 glaciers (covering 120 km^{2}) in our *V*–*A* relationship-based volume calculations.

## Data

### Calibration dataset for the volume–area relationship

To derive a reliable volume–area relationship specific to Svalbard glaciers, we need a large sample of accurately measured volume and area pairs, (*V*, *A*) pairs, representative of the size distribution of the total population of Svalbard glaciers. The individual glacier volumes are calculated from GPR-retrieved ice thickness data. Therefore, as a first step we compiled, under the framework of the European Science Foundation-supported PolarCLIMATE–SvalGlac project, an inventory of radio-echo sounded glaciers of Svalbard. This inventory is available online at http://svalglac.eu/. There are a total of 314 entries in the inventory, corresponding to 154 different glaciers (many glaciers have been surveyed more than once or have been surveyed using distinct radar equipment). Of these glaciers, we selected only 60 for our sample of (*V, A*) pairs. The selected sample consists of glaciers for which the net of GPR profiles covers most of the glacier basin and is dense enough to allow for a sufficiently accurate volume estimate.

Our sample of 60 (*V, A*) pairs was formed as follows. We calculated, from the original GPR ice thickness data, 36 glacier volumes with accuracy in general better than 10%. Of these 36 glaciers, 25 were radio-echo sounded during 1999–2014. Ten of them, located in western Nordenskiöld Land, are reported in Reference Martín-EspañolMartín-Español and others (2013), and another eight, located in Wedel Jarslberg Land (Fig. 1), are reported in Reference NavarroNavarro and others (2014). The remaining seven have not been published elsewhere; five of them are located in western and central Nordenskiöld Land (Fig. 1) and were echo sounded in spring 2013, while two are located in Sabine Land (Fig. 1) and were echo sounded in spring 2014. We also gathered from the literature the ice thickness maps and volumes of six other glaciers with a dense net of GPR profiles and accuracy better than 20% (according to the original sources). This totals 42 glaciers. Our procedure for estimating the error in glacier volume takes into account: (1) the point-dependent ice thickness errors of the GPR data (including the GPS-positioning error, the error in timing and the error in radio-wave velocity), (2) the interpolation error at every gridcell and (3) the volume error stemming from the uncertainty in the glacier boundary delineation. We estimate the interpolation error by using a modified ordinary kriging routine whereby the uncertainty in thickness at a given gridpoint is calculated relating the cross-validation errors with the distance to the nearest GPR measurement. We propagate the data errors to the gridpoints using the same weights adopted in the kriging interpolation procedure. Finally, the errors obtained at each gridpoint are optimally combined and averaged, considering the autocorrelation length scale of the ice thickness, to obtain the volume error. Further details are provided by Reference Martín-EspañolMartín-Español (2013).

A recent study (Reference Farinotti and HussFarinotti and Huss, 2013) has shown that, in estimating the accuracy with which the total volume of a glacier population can be recovered from a *V*–*A* relationship, the volume measurement uncertainty plays a secondary role compared with the sizes of both the total population and the sample, provided that the latter is sufficiently large. Following these ideas, we added to our sample some additional entries extracted from the catalogue of glacier volumes of the whole world (Reference Cogley, Henderson-Sellers and McGuffieCogley, 2012), reported in Reference GrinstedGrinsted (2013) (henceforth, referred to as Cogley’s catalogue). Within this catalogue, there are 31 glaciers in Svalbard and 19 of these are different from the 42 so far included in our set of (*V, A*) pairs. These 19 glaciers correspond to Scott Polar Research Institute–Norsk Polarinstitutt (Reference Dowdeswell, Drewry, Liestøl and OrheimDowdeswell and others, 1984) and Soviet (Reference YuYa and ZhuravlevMacheret and Zhuravlev, 1982) airborne radio-echo soundings in the 1970s and 1980s, which only covered the glacier centre lines, but exclude the glaciers with doubtful bed reflection interpretation from Soviet flights discussed by Reference Dowdeswell, Drewry, Liestøl and OrheimDowdeswell and others (1984) and later acknowledged by Reference Macheret, Zhuravlev and BobrovaMacheret and others (1984). These glaciers have been used by other authors (e.g. Reference GrinstedGrinsted, 2013) to build *V*–*A* relationships. Of these 19 glaciers we excluded one, Åsgårdfonna, because it is an ensemble of different glacier basins and it is not possible to calculate the volumes of its individual basins. With the addition of these 18 glaciers, for which we assume a volume accurate to better than 30% (Reference Martín-EspañolMartín-Español, 2013), we obtain our final number of 60 (*V, A*) pairs.

The basic data for the glaciers in our sample of 60 (*V, A*) pairs are given in the Appendix, and the location of the glaciers is shown in Figure 1.

### Comparison of the sample with the total population of Svalbard glaciers

The size distribution of our sample dataset ideally should be representative of the total population of Svalbard glaciers. The most up-to-date inventory of Svalbard glaciers is the Randolph Glacier Inventory V3.2 (Reference PfefferPfeffer and others, 2014), which contains the outlines for 1615 individual glacier basins in Svalbard (including Jan Mayen). According to this database, the total glacierized area of Svalbard is 33 922 km^{2}. This area reduces to 23 325 km^{2} if we exclude Austfonna and Vestfonna, Nordaustlandet (henceforth referred to as *NAL*).

Figure 2 shows the area distribution of the calibration dataset and compares it with the area distribution of the population of 1562 Svalbard and Jan Mayen glaciers excluding the glacier basins of Nordaustlandet’s ice caps. The RGI contains *∼*300 glaciers (20% of the inventory excluding *NAL*) larger than 10 km^{2}, which comprise 90% of the glacier area. The remaining 1261 glaciers (80% of the inventory excluding *NAL*) are <10 km^{2} and comprise 10% of the total glacier area. Among the latter, there are 610 glaciers, glacierets and snowpatches smaller than 1 km^{2}, representing 39% of the inventory and 1% of the total glacier area. In comparison, our sample of 60 glaciers has 31 glaciers (52% of the sample) larger than 10 km^{2}, comprising 95% of the total sample area, and only 29 glaciers (48% of the sample) equal to or smaller than 10 km^{2}, comprising 5% of the total area of the sample. Consequently, our sample dataset is biased towards large glaciers. This bias was anticipated, as small valley glaciers are customarily ignored when planning fieldwork because they are often difficult to access and seem less appealing than larger glaciers. To counterbalance this bias towards large glaciers in our sample dataset, we will later introduce an area-related weighting function.

Our calibration dataset is a fair representation of the total population of Svalbard glaciers in terms of terminus type (land-terminating or tidewater): 12 out of 60 glaciers in the sample (20%) are of tidewater type, while 9% of Svalbard glaciers are known to be of tidewater type, but this percentage increases to 14% when glaciers smaller than 1 km^{2} are excluded.

## Methods

We aim to calculate the total volume of Svalbard glaciers using scaling relationships. Volume–area scaling (e.g. Reference Chen and OhmuraChen and Ohmura, 1990; Reference Bahr, Meier and PeckhamBahr and others, 1997) allows estimation of the volume of a glacier from its area by means of a power law of the type

where *c* and *γ* are scaling parameters to be calibrated against a given set of (*V*, *A*) pairs, which can be done using different techniques, as described below.

We note that the value of *γ* derived by Reference Bahr, Meier and PeckhamBahr and others (1997), based on dimensional analysis, is *γ* = 1.375. Allowing slightly different *γ* values based on fits to observed data implies moving from a theoretical towards a statistical approach. This will become clearer later, when we introduce the multivariate scaling relations. Nevertheless, we remark that the *γ* value derived by Reference Bahr, Meier and PeckhamBahr and others (1997) from dimensional analysis is not entirely based on theoretical considerations, since it involves closure choices based on observational data. Moreover, the value *γ* = 1.375 is based on the assumption of a shallow-ice approximation dynamical model. Considering models slightly differing from it supports the possibility of slightly different values for the exponent *γ* in the *V*–*A* scaling relationship

### Regression techniques

We use two different regression techniques to estimate the parameters *c* and *γ* of the scaling relationships: (1) Least-squares regression in a log-log space (*logmse*) (Reference BahrBahr, 1997a):

where *n* is the total number of glaciers,*V*
_{model} are the volumes predicted by the scaling law with a set of parameters *p*, and *V*
_{obs} are the observed volumes in the glacier volume database. This model is very sensitive to outliers (see, e.g., Reference GrinstedGrinsted, 2013). (2) Least absolute deviation regression (*absdev*), proposed by Reference GrinstedGrinsted (2013), which minimizes the misfit function

where *A*
_{obs} is the observed area of each individual glacier. As noted by Reference GrinstedGrinsted (2013), this is a strategy best suited for sea-level rise studies, as it minimizes the absolute volume misfit, in addition to being robust to outliers and asymmetric distributions (Reference Cade and RichardsCade and Richards, 1996). This misfit function is weighted by the inverse of the square root of the area, which reduces the sampling bias towards larger glaciers previously identified in our calibration dataset.

### Scaling-law parameters for different glacier settings

It might be thought a priori that scaling relationships derived for specific types of glaciers would generate better volume predictions than a general-purpose *V*–*A* scaling relationship. However, their expected accuracy improves both with the size of the total target population of glaciers and with the size of the sample used to derive the scaling-law parameters (Reference Farinotti and HussFarinotti and Huss, 2013). Partitioning a given sample into subgroups by specific characteristics such as glacier size, shape or slope reduces the size of the sample from which the scaling parameters are derived, and hence reduces the expected accuracy in the total volume estimate. Therefore there is a trade-off between the improvement expected by the use of glacier-type specific scaling relationships and the worsening yielded by the reduction of the sample size. We undertook an experiment aimed at verifying whether scaling relationships obtained through characterization of individual glacier size or morphology imply significant differences in the estimated volume of Svalbard glaciers, and whether this partitioning implies any noticeable pattern in the scaling exponents, indicative of the influence of particular glacier settings. In particular, we considered partitions of our available sample of Svalbard glaciers by size, shape and slope, as done by Reference Adhikari and MarshallAdhikari and Marshall (2012). Shape refers to the horizontal characteristic glacier shape, given by *W/L*, where *L* is the glacier length measured along the central flowline and *W* is the mean glacier width, defined as *W = AIL*, where *A* is glacier area. Thus our measure of shape is *A/L*
^{2}. Slope is the mean bedrock slope in the principal flow direction. As it is not available from measurements done at the glacier surface, following Reference Adhikari and MarshallAdhikari and Marshall (2012) we took the mean surface slope, given by *E/L*, where *E* is the elevation range, as a proxy of the mean bedrock slope.

### Multivariate analysis

Multivariate analysis has proved successful for strengthening the capacity of scaling relationships to estimate glacier volume (Reference GrinstedGrinsted, 2013). We adopt a multivariate approach to predict the total volume of Svalbard glaciers from a combination of different predictors: glacier area (*A*), glacier maximum length (*L*) and glacier elevation range (*E*). We excluded from this analysis variables such as glacier slope and shape to avoid multicollinearity. We note that, by including additional predictors, we leave the realm of physically based scaling and move towards a statistically based relation. Thus, from statistical tests pursuing the analysis of the variance we found the most significant variables to which to fit the distinct scaling models.

### Error estimates

The volume estimates based on *V*–*A* scaling relationships can involve very large errors when applied to individual glaciers, though the errors are much lower, because of statistical compensation, when applied to large sets of glaciers. According to Reference MeierMeier and others (2007), estimated volume errors for individual glaciers could exceed 50% but these uncertainties are reduced to 25% for an ensemble of glaciers. Consistent with these results, a recent study by Reference Adhikari and MarshallAdhikari and Marshall (2012), using a *V*–*A* relationship based on a sample of 280 synthetic random mountain glaciers, has shown that, when estimating the volumes for all individual glaciers in their ensemble, the average glacier volume error is small (2.8%), with a mean absolute error of 18.3%. Their interquartile spread of results is also reasonable, with 50% of errors in the range [–14.7, 15.9]%, but errors for individual glaciers can be high, in the range [–47.6, +99.7]%. Their mean and mean absolute errors for individual glacier volume estimates are reduced to 1.5% and 14.4%, respectively, when shape-based parameters are applied to their random ensemble of glaciers. But these error estimates are too optimistic for a real case in a sense: the sizes of the total population of glaciers to which the scaling relationship is applied and of the calibration dataset from which the scaling parameters are derived are equal, which does not occur in real applications.

An alternative approach is that of Reference Farinotti and HussFarinotti and Huss (2013), who have derived upper bound estimates for the accuracy with which the total volume can be recovered in relation to the size of the glacier population, the size of the sample of glaciers used to estimate the scaling parameters and the uncertainty of the measured volumes. Their study has shown that: (1) a low level of accuracy is expected if scaling is applied with coefficients estimated from a small set of samples; (2) accuracy increases with increasing size of the considered glacier population; and (3) volume measurement uncertainty plays a secondary role if a sufficient number of measured glacier volumes are available. The last of these ideas was behind our decision to add to our initial calibration dataset of 42 glaciers, with volume accuracies <20%, 18 additional glaciers with volume accuracies <30%. According to the estimates of Reference Farinotti and HussFarinotti and Huss (2013), when computing the total volume of ∼1500 Svalbard glaciers using a *V*–*A* relationship based on our sample of 60 measured (*V*, *A*) pairs, whose measured volumes have typical accuracies within 5–25%, we could expect a priori an accuracy with an upper bound of 45–50% (Fig. 3). This is an upper bound for the accuracy, and thus a lower bound for the error, derived from synthetic data and only to be reached under ideal conditions (Reference Farinotti and HussFarinotti and Huss, 2013).

We note, however, that the study by Reference Farinotti and HussFarinotti and Huss (2013) is based on the assumption of independent and identically distributed (i.i.d.) data pairs. The individual volumes of the global population of glaciers, and the measured volumes and areas, are all assumed to have errors that follow normal (Gaussian) distributions with zero mean and given standard deviations. These assumptions are violated in our regional study, so there would be no contradiction if our Svalbard-specific scaling law, based on a rather large sample of measured glacier volumes, representative (in terms of spatial distribution and morphology) of the total population of Svalbard glaciers, produced better results than those predicted by Reference Farinotti and HussFarinotti and Huss (2013). We thus adopt an alternative approach, and estimate the error of our total volume calculation following the principles of error propagation (Reference Bevington and RobinsonBevington and Robinson, 2003) and the rationale of Reference Radić and HockRadić and Hock (2010). From the *V*–*A* scaling model, the total volume is computed as

The variables and parameters involved in the error in *V* are thus *Ai*, *i* 1/4 1, …, *n*, *c* and *γ*, and we assume that they are random and independent. Then, from error propagation,

where *δc*, *δγ* and *δA _{i}
* represent the errors in the corresponding parameters and variables. As typical error in area for Svalbard glaciers we take the value of 8% suggested by Reference König, Nuth, Kohler, Moholdt, Pettersen, Kargel, Leonard, Bishop, Kääb and RaupKönig and others (2014). Following the procedure suggested by Reference BahrBahr (1997b), we estimate the error in

*c*as follows. We set a fixed value of

*γ*(the value obtained for the

*V*–

*A*regression retrieved using the 60 glaciers in our sample) and compute 60 individual values for

*c*using . We then take as error in

*c*the standard deviation of the distribution of

*ck*(Fig. 4a), which is

*c*= 0

*:*0109. Likewise, for estimating the error in

*γ*we set a fixed value of

*c*(the value obtained for the

*V*–

*A*regression retrieved using the 60 glaciers in the sample) and compute 60 individual values for

*γ*using . We then take as error in

*γ*the standard deviation of the distribution of

*γ*(Fig. 4b), which is

_{k}*γ*= 0

*:*229. Using the equations and values above, together with the RGI V3.2 for Svalbard (excluding Nordaustlandet ice caps), gives relative errors in the total volume estimate of 767 and 891 km

^{3}, when the

*logmse*and

*absdev*fitting strategies, respectively, are used to derive the

*V*–

*A*relationships. The procedure for estimating the error in the scaling exponent

*γ*described above differs from the approach taken by Reference Radić and HockRadić and Hock (2010), who approximate

*δγ*by the difference between the value

*γ*= 1.375 derived theoretically by Reference Bahr, Meier and PeckhamBahr and others (1997) and that obtained from the

*V*–

*A*fit to the observations (Reference BahrBahr, 1997b). Our procedure gives, for the Svalbard case, a larger but more realistic error estimate of

*γ*(and hence of the volume), and is consistent with the procedure used for estimating the error in the scaling coefficient

*c*.

## Results and Discussion

### Regression techniques

The results of the calibration of the simple *V*–*A* model given by Eqn (1) using our dataset of 60 (*V*, *A*) pairs, by means of both *logmse* and *absdev* regression techniques, are shown in Table 2 and Figure 5. Table 2 also includes the cross-validation error based on the calibration dataset, estimated by leave-one-out cross-validation (Reference GeisserGeisser, 1993), and the coefficient of determination *R*
^{2} of the fit, and also presents the volumes computed by each scaling relationship and their estimated errors, calculated as described in the previous subsection. The absolute errors in volume shown in Table 2 correspond to relative errors of 20% (*logmse*) and 24% (*absdev*), and the volume estimates differ from each other by ∼5%.

Cross-validation is estimated as follows. In *k*-fold cross-validation, the original sample is randomly partitioned into *k* equal-size subsamples. Of the *k* subsamples, a single subsample is retained as the validation data for testing the model, and the remaining *k*-1 subsamples are used as training data. The cross-validation process is then repeated *k* times, with each of the *k* subsamples used exactly once as validation data. Leave-one-out cross-validation is the same as a *k*-fold cross-validation with *k* being equal to the number of observations in the original sample. In each case, all observations are used for both training and validation, and each observation is used for validation exactly once. The cross-validation errors shown in Table 2 are the root-meansquare errors (RMSE) for the volumes estimated during the leave-one-out cross-validation procedure.

To verify that our derived scaling relationships are not influenced by the original bias in the size distribution of the glacier sample, we performed, for each of the regression techniques, the following test. For each particular glacier in our sample, we computed the residual volume obtained by subtracting, from the volume of the glacier computed from GPR data, the estimated volume of that particular glacier derived from a *V*–*A* scaling law obtained by removing that particular glacier from the sample. We then plotted the residual volume versus the area of the glacier. No significant correlation was found, indicating that there is no noticeable influence of the original bias in the size distribution of the glacier sample on our derived scaling relationships.

### Scaling-law parameters for different glacier settings

The results for the best-fitting scaling parameters *c* and *γ* for specific subgroups of glaciers arranged by size, shape and slope, using both the *logmse* and *absdev* fitting strategies, are shown in Table 3. Overall, partitioning by size yields poorer fits and the smallest volume for the entire population of Svalbard glaciers excluding *NAL* (*V* ⊂ [3561, 3656] km^{3}), whereas the largest correspond to the partitionings by shape (*V* ⊂ [3867, 3897] km^{3}) and slope (*V* ⊂ [3579, 4144] km^{3}).

We also calculate *γ* for a constrained value of *c* = 0*:*0343, taken from the *logmse* regression for the entire population of glaciers (Table 2). Here we use the *logmse* misfit function rather than the *absdev* approach because *logmse* is most often used in the literature and provides a value for *c* closer to those obtained by other authors, which facilitates comparison of the *γ* values. The results of the constrained experiment show minimum *γ* values for the highest ranges of slope and shape, in agreement with the results obtained by Reference Adhikari and MarshallAdhikari and Marshall (2012) for a collection of synthetic glaciers. This means that the volumes of steep-slope and cirque-type glaciers are less sensitive to changes in glacier area. This could actually be an indication that response timescales play a role. In particular, long, steep glaciers presumably have a different response time than short, flat glaciers, which means that current disequilibrium between volume and shape, caused by the ongoing mass losses, expresses itself in different optimal parameter values. Recalling, from Table 2, that the total volume of Svalbard glaciers excluding *NAL*, computed using a single scaling law for all glaciers and the *logmse* misfit function, was 3955 km^{3}, we see that the constrained experiment gives volumes that are ∼2% smaller when using partitionings by size or shape, while nearly equal (only slightly larger) when using partitioning by slope. We note, however, the small number of glaciers belonging to some of the subgroups, which probably lie near the minimum admissible number for reaching statistically significant conclusions.

In our study of Svalbard glaciers, a straightforward partitioning would be to distinguish between tidewater and land-terminating glaciers. However, for our sample of (*V*, *A*) pairs this partitioning is nearly equivalent to a partitioning by size (with the tidewater glaciers being the largest), so it gives no further insight. Another clear classification criterion would be a grouping into surging and non-surging glaciers. In this case, the available set of (*V*, *A*) pairs for surging glaciers is so small that it does not allow us to derive a specific scaling law. Moreover, the volume–area ratios of surging glaciers are dependent upon the surge phase: while during the early post-surge period surging glaciers will generally have smaller volume–area ratios than non-surging glaciers, during the build-up period of the surge the converse will generally be true. As a sensitivity test, we removed the surge glaciers from our calibration dataset and, keeping *c* fixed, we recalculated the value of the exponent *γ*. This resulted in a slightly larger scaling exponent, *γ* = 1.333 (using the *logmse* fitting strategy), which indicates that most of the surge glaciers included in our sample are in their early post-surge period, with a consequently smaller volume–area ratio.

### Multivariate analysis

A multivariate statistical analysis reveals that a simple model *V* ∝ *A ^{γ}
* in the log-log space explains 98.6% of the variance. However, the analysis of the cross-validation errors reveals that results are improved when all the variables are included in the model, i.e. , especially for the

*logmse*misfit function, as more significant glaciological information is used to predict glacier volume. Table 4 shows the ice volume estimates for Svalbard using different regression models, and their accuracy, estimated in terms of leave-oneout cross-validation errors. Accuracy is strongly reduced in the models where the variable

*A*is not present (14.9% and 22.9% cross-validation RMSE), in contrast to the results obtained by Reference Radić, Hock and OerlemansRadić and others (2008), where

*V*∝

*L*had the lowest error. In general the

^{γ}*logmse*strategy provides lower cross-validation errors, the lowest corresponding to the scaling law, although the small exponent of

*E*(Table 4) confirms the low predicting capability of this variable.

### Comparison with other *V*–*A* scaling relationships

In Table 5 we compare our results for the glacier volume of Svalbard excluding Nordaustlandet (*V*
_{Svexcl:NAL}) with those obtained using other scaling relationships available in the literature. For the sake of homogeneity, the RGI V3.2 is used for all of them, and we use our standard *V*–*A* scaling law (instead of the multivariate), as most published relationships are of the *V*–*A* type. Also for this reason we choose, among the distinct relationships presented by Reference GrinstedGrinsted (2013) (several of them multivariate), his standard *V*–*A* scaling law. All scaling relationships included in the table are of the *V*–*A* type, except those of Reference Hagen, Liestøl, Roland and JørgensenHagen and others (1993) and Reference Huss and FarinottiHuss and Farinotti (2012), which are of the *H*–*A* type, with *H* the mean glacier thickness. Because of this, the exponent *γ* in Huss and Farinotti’s relationship is nearly 1 unit lower than those for the univariate *V*–*A* relationships, and also the coefficient *c* is substantially different. As a complementary test of performance, we also present the percentage of error incurred when calculating the total volume of our calibration dataset using the relationships from other published scaling models (Table 5; we exclude ours, since their results would be skewed). The estimates in Table 5 range within 2462–6170 km^{3}, with an average value of 4421 ± 1037 km^{3} if we exclude our own estimates, and 4334 ± 971 km^{3} if our results are included. In both cases, the quoted errors indicate the standard deviations of the different estimates considered, which illustrates the spread of the results.

In addition to the scaling relationships derived in this study, those by Reference Macheret, Zhuravlev and BobrovaMacheret and others (1984) and Reference Hagen, Liestøl, Roland and JørgensenHagen and others (1993) are specific to Svalbard glaciers. However, those two studies were based on a glacier inventory of Svalbard several decades older than ours, so their areas and volumes are generally larger due to the overall retreating and thinning trends of Svalbard glaciers during recent decades (Reference Nuth, Kohler, Aas, Brandt and HagenNuth and others, 2007, Reference Nuth, Moholdt, Kohler, Hagen and Kääb2010). As shown in Table 5, the scaling approach of Reference Macheret, Zhuravlev and BobrovaMacheret and others (1984) overestimates the total volume of the calibration dataset by 32%, while that of Reference Hagen, Liestøl, Roland and JørgensenHagen and others (1993) also overestimates it, but only by 9%. This could be an indication that the volume–area relationship of Svalbard glaciers has varied over the past 30 years, with mass loss dominated by glacier thinning rather than front retreat, which is consistent with observations (e.g. Reference Nuth, Kohler, Aas, Brandt and HagenNuth and others, 2007, Reference Nuth, Moholdt, Kohler, Hagen and Kääb2010), resulting in smaller *V/A* ratios. But these results could also partly reflect some bias in the glacier volume data used to derive the mentioned scaling relationships.

By far the closest estimates to both our results for the total volume of Svalbard glaciers (excluding *NAL*) and the total volume of the calibration dataset are those obtained using the scaling relationships of Reference Chen and OhmuraChen and Ohmura (1990) and Reference BahrBahr (1997a). These scaling approaches are based on a global sample of valley glaciers. We might think that, having excluded Austfonna and Vestfonna, the calibration datasets of both these global relationships and our regionally based one are somehow similar. However, the good fit is still surprising considering the large proportion of tidewater glaciers within the sample of Svalbard glaciers compared with the samples by Reference Chen and OhmuraChen and Ohmura (1990) and Reference BahrBahr (1997b). The scaling approaches of Reference Huss and FarinottiHuss and Farinotti (2012) and Reference GrinstedGrinsted (2013) also provide total volume estimates close to our own. However, when estimating the volume of the calibration dataset, the relationship by Reference Huss and FarinottiHuss and Farinotti (2012), which is a regional scaling law derived from the entire set of Svalbard glacier volumes computed with the physically based method developed by those authors, gives a large overestimate of 47%. A possible explanation is that, being a scaling relation based on the complete set of volumes, it gives a good approach to the volume of the total population, while, when applied to our sample, it does not provide such a good fit because the sample is not fully representative of the size distribution of the complete population of Svalbard glaciers.

The largest estimates of the total volume of Svalbard excluding *NAL* are those obtained using the scaling by Reference Adhikari and MarshallAdhikari and Marshall (2012), derived from their entire set of synthetic glaciers (without any partitioning by glacier types), and that based on partitioning by size. These also produce the largest overestimates of the volume of the calibration dataset. By contrast, the relationship based on partitioning by slope gives by far the lowest total volume estimate (but overestimates the volume of the sample), while that based on partitioning by shape produces a mid- to high total volume estimate (but underestimates the volume of the sample). Table 5 shows that, in general, a high total volume estimate is accompanied by an overestimate of the volume of the sample, and conversely. The scalings of Reference Adhikari and MarshallAdhikari and Marshall (2012) based on partitionings by slope and shape are an exception to this rule. The most likely reason is the irregular distribution of glaciers in our sample among the three slope and shape ranges considered by Reference Adhikari and MarshallAdhikari and Marshall (2012), 50–10–0 (slope) and 16–35–9 (shape), while the distribution is more regular in the case of partitioning by size (26–10–24). We note that the ranges of the partitionings by size, shape and slope in our experiment in the subsection ‘Scaling-law parameters for different glacier settings’ above, and in Reference Adhikari and MarshallAdhikari and Marshall (2012) are different. This is why our experiment using partitioning (which had a fair distribution of glaciers among the different ranges of variation of a given attribute) produced good results (consistent with the experiments without partitioning), while Reference Adhikari and MarshallAdhikari and Marshall (2012) also produced good results when working with their sample of glaciers (which also had a fair distribution of glaciers among the different ranges of each partitioning), but not when applied to our population and sample of Svalbard glaciers. From this we conclude that the scaling relationships based on partitionings by glacier attributes such as size, shape or slope are expected to produce good results provided the size of the sample is large enough and has a fair distribution of glaciers among the different ranges of values of each attribute, which should also be representative of the distribution existing in the total glacier population. While this representativeness calls for regionally based relationships, the main difficulty faced in this case is the requirement of a large sample. Our Svalbard partitioning experiment seems to be near the minimum admissible number of sampled glaciers.

The total volume estimate obtained using the scaling by Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others (2014) is also large, and their relationship greatly overestimates (by 43%) the volume of the calibration dataset. This illustrates how sensitive are the volume calculations from volume–area relationships to the choice of the scaling parameters, since Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others (2014) adopt a scaling law that uses the *c* coefficient obtained by Reference Chen and OhmuraChen and Ohmura (1990) (which produced results close to ours when used with its own exponent *γ* = 1*:*357) together with the exponent *γ* = 1*:*375 derived by Reference Bahr, Meier and PeckhamBahr and others (1997) based on physical considerations for mountain glaciers (which we note is different from the *γ* = 1*:*36 obtained by Reference BahrBahr (1997a) that gave results close to ours when used with its associated *c* coefficient). It should be noted that the *c* coefficients in the *V* –*A* scaling laws in Table 5 are given in units of km^{3−2γ
} and therefore depend on the choice of *γ*. Since the fits for *c* and *γ* are obtained simultaneously, and the parameter values are therefore mutually influenced, we do not recommend taking separate values for *c* and *γ* from independent studies. If a *γ* value is taken from a theoretical study, as done by Reference Radić, Bliss, Beedlow, Hock, Miles and CogleyRadić and others (2014) following Reference Radić and HockRadić and Hock (2010), then the *c* coefficient should ideally be obtained from a fit to measured volumes (keeping fixed the selected *γ* exponent).

### Total volume of Svalbard glaciers and potential contribution to sea-level rise

Our best estimates of the total volume of Svalbard glaciers, excluding Austfonna and Vestfonna, are those given by our multivariate scaling relationships including the glacier area, length and altitude range as variables. The results obtained using the *logmse* and *absdev* regression techniques, extracted from Table 4, are given in Table 6 accompanied by their error estimates calculated as described in the ‘Error estimates’ subsection above. To obtain the total volume of Svalbard glaciers we simply add the volumes of Austfonna and Vestfonna determined from extensive radio-echo sounding (both airborne and surface-based) as described by Reference DunseDunse (2011), Reference Pettersson, Christoffersen, Dowdeswell, Pohjola, Hubbard and StrozziPettersson and others (2011) and Reference Martín-EspañolMartín-Español (2013), which are 2559 km^{3} (Austfonna) and 442 km^{3} (Vestfonna), with respective relative errors of 3% and 7%, estimated by Reference Martín-EspañolMartín-Español (2013). Thus, the ice volume of Nordaustlandet ice caps totals 3001 ± 81 km^{3}. The resulting total volume of Svalbard glaciers is given in Table 6, where the errors from scaling (Svalbard excluding *NAL*) and GPR (*NAL*) are combined as the square root of the sum of squares. The values based on both regression techniques are shown, and we can consider their average, 6700 ± 835 km^{3}, as our best estimate of Svalbard ice volume. In terms of sea-level equivalent (SLE), assuming an oceanic area of 3.62 × 10^{8} km^{2} and a glacier ice density of 900 kg m^{−3}, our volume estimate corresponds to a total potential contribution to sea-level rise of 17 ± 2 mm SLE.

This volume estimate is in the low range of those published in the literature (Table 1), including the regionally based estimates of Reference Macheret, Zhuravlev and BobrovaMacheret and others (1984) and Reference Hagen, Liestøl, Roland and JørgensenHagen and others (1993), but higher than that of Reference GrinstedGrinsted (2013). We note that the volumes given in Table 1 are based on different inventories. This should have an impact especially on the estimates by Reference Macheret, Zhuravlev and BobrovaMacheret and others (1984) and Reference Hagen, Liestøl, Roland and JørgensenHagen and others (1993). The most recent estimates (including ours) are all based on the RGI. Although the RGI version used sometimes differs, this should have little impact on the results. Although the number of glacier complexes worldwide is different in RGI V1.0 and V2.0 compared with V3.2 (Reference PfefferPfeffer and others, 2014), and this could have an impact on the volume calculations, the number of glacier complexes in the particular case of Svalbard has not changed between RGI versions (Reference ArendtArendt and others, 2014). The only change from V1.0 to V2.0 that has an impact on glacier volume is that the outlines of Jan Mayen were included. This implies an increase in volume of ∼8–9 km^{3} (depending on the scaling relation used), which only applies to the comparison between the estimate by Reference Marzeion, Jarosch and HoferMarzeion and others (2012), based on RGI V1.0, and all later ones, based on V2.0 except ours, based on V3.2. We also recall that all results except that of Reference Huss and FarinottiHuss and Farinotti (2012) (derived from a physically based method) are based on scaling relationships, and that our estimate combines scaling (for Svalbard excluding Nordaustlandet) and direct calculation from GPR-retrieved ice-thickness data (for Nordaustlandet ice caps). We also note that, in all versions of the RGI, Austfonna and Vestfonna (and also ice caps covering other islands of the Svalbard archipelago) appear subdivided into individual drainage basins, and that currently the *RGIflag* attribute of the RGI only distinguishes between glaciers and ice caps for the Antarctic and Subantarctic RGI region (Reference PfefferPfeffer and others, 2014). Consequently, any scaling-based estimate, even if it has separate scaling laws for glaciers and ice caps, if applied to Svalbard (or any region other than the Antarctic and Subantarctic) using the RGI without corrections, will compute the volume of each entry in the inventory as if it were a glacier. This will likely lead to biased results. With the aim of quantifying the differences for the case of Austfonna and Vestfonna, we calculated their volumes as the sum of the volumes of all their basins, as defined in the RGI V3.2, using our *V*-*A* scaling (as done, e.g., by Reference GrinstedGrinsted, 2013). This gives a total volume of Nordaustlandet ice caps of 2444 km^{3} (averaging the *logmse*- and *absdev*-based results), which is 557 km^{3} lower than the volume of 3001 km^{3} obtained from the GPR data. Nearly all the difference (546 km^{3}) corresponds to Austfonna. Using the scaling-based estimate for all of Svalbard gives a total volume of 6143 km^{3}, making our result closer to that of Reference GrinstedGrinsted (2013). It would be even closer (6003 km^{3}) if we took our *absdev*-based result (recall that Grinsted’s result is also based on the *absdev* regression). We stress, anyway, that the GPR-based volume estimate for Austfonna and Vestfonna is more credible than any estimate based on scaling relationships.

## Summarizing Conclusions

Using 60 highly accurate Svalbard volume–area pairs, with volumes calculated from GPR data, we calibrated different scaling relationships with the aim of estimating the total volume of Svalbard glaciers and, thus, their potential contribution to sea-level rise. These relationships were applied to each glacier record included in the RGI V3.2 for Svalbard and Jan Mayen. We estimated the total ice volume of Svalbard glaciers at 6700 ± 835 km^{3}, or 17 ± 2 mm SLE. We obtained evidence, from the slope- and shape-based scaling relationships, that the volumes of steep-slope and cirque-type glaciers appear to be less sensitive to changes in glacier area. An alternative interpretation is that they could just have different response times. The multivariate analysis shows minimum cross-validation errors when the volume– area–length–elevation range scaling model is used, for both the *logmse* and *absdev* regression techniques, which suggests that the use of multivariate scaling relationships improves the accuracy as compared with the standard volume–area scalings. Our estimate of the Svalbard glacier volume lies in the low range of previously published estimates, with only the estimate by Reference GrinstedGrinsted (2013) providing a lower volume. However, a fair comparison is difficult, as the published estimates are based on different inventories. Fortunately, the most recent estimates are all based on the RGI, and, even if based on different versions of it, the only difference corresponds to the 8–9 km^{3} of added volume corresponding to the Jan Mayen glaciers from RGI V2.0 onwards. An additional difficulty is that not all scaling relations are of the *V* –*A* type, but some include further variables. Also, the way of dealing with Austfonna and Vestfonna ice caps, which appear in the RGI subdivided into individual basins, all classified as *glaciers*, could make a difference. For this reason, we made a comparison based on the use of the RGI V3.2 for all scaling relationships, and restricted to Svalbard glaciers excluding Austfonna and Vestfonna. This comparison showed that globally based relationships sometimes provide results close to regionally based ones. Their suitability, rather than depending on being globally or regionally based, depends on how well the sample of glaciers from which the scaling relationship is derived represents the actual distribution of the total population of glaciers of the region under study. Of course other factors, such as the sample size, the total population size and the accuracy of the individual volumes in the sample, also play an important role, as shown by Reference Farinotti and HussFarinotti and Huss (2013). Regarding this, our calibration dataset of 60 (*V, A*) pairs, with volume errors typically within 5–25% but quite often <10%, and being a fair representation of the size distribution of the total population of Svalbard glaciers, strongly supports the accuracy of our estimate in comparison with other published estimates. The accuracy of our estimate is also expected to have been improved by the use, for Austfonna and Vestfonna, of the volume calculated directly from the GPR-retrieved ice-thickness data. The comparison also demonstrated the sensitivity of the volume estimates to the scaling parameters and, in particular, that the coefficient *c* and the exponent *γ* should be fitted simultaneously. More specifically, if one of these is taken from theoretical considerations, the other should be fitted against available volume measurements rather than taking it from other fits available in the literature. Finally, a note of caution is sounded about the use of the RGI together with scaling relationships distinguishing between glaciers and ice caps, since none of the currently available RGI versions identifies the ice caps (i.e. all inventory entries appear as *glacier*), except for the Antarctic and Subantartic RGI regions. This calls for the completion of this important attribute in the RGI.

## Acknowledgements

This research was supported by grant EUI2009-04096 (PolarCLIMATE-SvalGlac) from the Spanish EuroResearch Programme, grants CGL2005-05483, CTM2008-05878 and CTM2011-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 grants 10-05-00133-a and 11-05-00728-a from the Russian Fund of Basic Research. We are grateful for the suggestions by Graham Cogley, Surendra Adhikari, Ben Marzeion and an anonymous reviewer, as well as the work of the scientific editor, Jo Jacka, which greatly improved the manuscript.

## Appendix