Estimate of the total volume of Svalbard glaciers, and their potential contribution to sea-level rise, using new regionally based scaling relationships

Abstract

We present a set of new volume scaling relationships specific to Svalbard glaciers, derived from a sample of 60 volume–area pairs. Glacier volumes are computed from ground-penetrating radar (GPR)-retrieved ice thickness measurements, which have been compiled from different sources for this study. The most precise scaling models, in terms of lowest cross-validation errors, are obtained using a multivariate approach where, in addition to glacier area, glacier length and elevation range are also used as predictors. Using this multivariate scaling approach, together with the Randolph Glacier Inventory V3.2 for Svalbard and Jan Mayen, we obtain a regional volume estimate of 6700 ± 835 km³, or 17 ± 2 mm of sea-level equivalent (SLE). This result lies in the mid- to low range of recently published estimates, which show values as varied as 13 and 24 mm SLE. We assess the sensitivity of the scaling exponents to glacier characteristics such as size, aspect ratio and average slope, and find that the volume of steep-slope and cirque-type glaciers is not very sensitive to changes in glacier area.

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 (Stocker 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% (Radić 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 (Radić and Hock, 2010; Huss and Farinotti, 2012; Marzeion and others, 2012; Grinsted, 2013; Radić 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 Cogley (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 (Pfeffer 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 (Bahr 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. McNabb 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 (Bahr 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 (Radić and Hock, 2010; Marzeion and others, 2012; Grinsted, 2013; Radić 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 (McNabb and others, 2012)), physically based methods are often restricted to the analysis of unique or small sets of glaciers (e.g. Farinotti and others, 2009). However, recent applications have dealt with much larger sets of glaciers, at both regional (Frey and others, 2013; Huss and Farinotti, 2014) and global (Huss 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 (Huss 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 Radić 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.

Fig. 1. Radio-echo sounded glaciers on Svalbard for which a GPR-based volume estimate is available. Each glacier is numbered as in the Appendix, using the IDs taken from the inventory of radio-echo sounded glaciers available online at http://svalglac.eu/. The basic data for each glacier are also given in the Appendix. Austfonna and Vestfonna ice caps appear subdivided into their individual basins, as given in the RGI, although they are not used in our V–A relationships. The numbers for the glaciers echo sounded by our research team and/or colleagues appear with white background, and all have estimated relative errors in volume <10%. The rest of the glaciers have been echo sounded by other authors. Those with orange background have relative errors in volume <20% (about half of them <10%), while we assume that those with yellow background, taken from Cogley’s catalogue as given in Grinsted (2013), have estimated errors in volume <30%.

The ice masses of Svalbard cover ∼33 922 km² (Pfeffer 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 (Marzeion and others, 2012; Grinsted, 2013; Radić and others, 2014) and physically based methods as in Huss and Farinotti (2012), revealing high variance across different estimates (Table 1). Older estimates derived from Svalbard-specific scaling relationships (Macheret and Zhuravlev, 1982; Hagen 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 (Macheret and Zhuravlev, 1982; Dowdeswell 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.

Table 1. Total volume estimates for Svalbard glaciers (V _Sv) found in the literature. All estimates from 2012 onwards use the Randolph Glacier Inventory V2.0, except Marzeion and others (2012), which uses V1.0

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 (Pfeffer 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²) 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²) 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 Martín-Español and others (2013), and another eight, located in Wedel Jarslberg Land (Fig. 1), are reported in Navarro 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 Martín-Español (2013).

A recent study (Farinotti 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 (Cogley, 2012), reported in Grinsted (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 (Dowdeswell and others, 1984) and Soviet (Macheret 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 Dowdeswell and others (1984) and later acknowledged by Macheret and others (1984). These glaciers have been used by other authors (e.g. Grinsted, 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% (Martí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 (Pfeffer 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². This area reduces to 23 325 km² 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², which comprise 90% of the glacier area. The remaining 1261 glaciers (80% of the inventory excluding NAL) are <10 km² and comprise 10% of the total glacier area. Among the latter, there are 610 glaciers, glacierets and snowpatches smaller than 1 km², 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², comprising 95% of the total sample area, and only 29 glaciers (48% of the sample) equal to or smaller than 10 km², 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.

Fig. 2. Area distribution of the sample dataset versus area distribution of the complete population of Svalbard glaciers, excluding NAL.

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² are excluded.

Methods

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

(1)

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 Bahr 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 Bahr 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) (Bahr, 1997a):

(2)

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., Grinsted, 2013). (2) Least absolute deviation regression (absdev), proposed by Grinsted (2013), which minimizes the misfit function

(3)

where A _obs is the observed area of each individual glacier. As noted by Grinsted (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 (Cade 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 (Farinotti 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 Adhikari 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 ². Slope is the mean bedrock slope in the principal flow direction. As it is not available from measurements done at the glacier surface, following Adhikari 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 (Grinsted, 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 Meier 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 Adhikari 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 Farinotti 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 Farinotti 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 (Farinotti and Huss, 2013).

Fig. 3. Expected accuracy with which the total true volume of a population of 1500 synthetic glaciers will be recovered (at 95% confidence) using volume–area scaling relationships. The different curves correspond to different measurement uncertainties of individual glacier volumes (UKV). This figure was kindly generated by D. Farinotti.

We note, however, that the study by Farinotti 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 Farinotti and Huss (2013). We thus adopt an alternative approach, and estimate the error of our total volume calculation following the principles of error propagation (Bevington and Robinson, 2003) and the rationale of Radić and Hock (2010). From the V–A scaling model, the total volume is computed as

(4)

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,

(5)

(6)

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 König and others (2014). Following the procedure suggested by Bahr (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 γ_k (Fig. 4b), which is γ = 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³, 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 Radić and Hock (2010), who approximate δγ by the difference between the value γ = 1.375 derived theoretically by Bahr and others (1997) and that obtained from the V–A fit to the observations (Bahr, 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.

Fig. 4. Histograms of the distributions for parameters γ and 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 (Geisser, 1993), and the coefficient of determination R ² 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%.

Fig. 5. Volume–area scaling models calibrated using the logmse and absdev misfit functions with the 60 (V, A) pairs of Svalbard glaciers. The parameter set is given in Table 2.

Table 2. Derived scaling law and associated total volume of Svalbard glaciers (excluding NAL) calculated using the logmse and absdev regression techniques. The cross-validation error ∊ _crossval and the coefficient of determination R ² are also given.

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³), whereas the largest correspond to the partitionings by shape (V ⊂ [3867, 3897] km³) and slope (V ⊂ [3579, 4144] km³).

Table 3. Scaling-law parameters for different glaciological settings (following Adhikari and Marshall (2012), adapting ranges to our sample distribution). In addition to the logmse and absdev fitting strategies, results are also presented for a constrained experiment in which a fixed value of c is used (c = 0:0343) and only γ is fit. The coefficients of determination R ² of the different regressions are given, as well as the total volume of Svalbard glaciers V _Sv (excluding NAL) calculated using the given partitions into subgroups of glaciers. Values in brackets indicate the number of glaciers in each subgroup

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 Adhikari 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³, 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 Radić 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.

Table 4. Scaling laws resulting from the multivariate analysis, and their associated estimates of the total volume of Svalbard (excluding NAL), together with the cross-validation errors incurred when calculating the total volume of the 60 glaciers in our calibration dataset using these scaling laws, and the corresponding coefficients of determination R ²

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 Grinsted (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 Hagen and others (1993) and Huss 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³, with an average value of 4421 ± 1037 km³ if we exclude our own estimates, and 4334 ± 971 km³ 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.

Table 5. Total volume estimates for Svalbard glaciers (V _Sv), excluding NAL, calculated using various scaling relationships found in the literature, together with our own, and relative errors produced when these scaling laws are applied to our calibration dataset (∊V _Sample). Note that the values of the coefficient c in the V–A scaling laws are given in km^3−2γ , with γ the exponent of the scaling law, while in the literature they are often given in m^3–27γ

In addition to the scaling relationships derived in this study, those by Macheret and others (1984) and Hagen 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 (Nuth and others, 2007, 2010). As shown in Table 5, the scaling approach of Macheret and others (1984) overestimates the total volume of the calibration dataset by 32%, while that of Hagen 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. Nuth and others, 2007, 2010), 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 Chen and Ohmura (1990) and Bahr (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 Chen and Ohmura (1990) and Bahr (1997b). The scaling approaches of Huss and Farinotti (2012) and Grinsted (2013) also provide total volume estimates close to our own. However, when estimating the volume of the calibration dataset, the relationship by Huss 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 Adhikari 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 Adhikari 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 Adhikari 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 Adhikari 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 Adhikari 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 Radić 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 Radić and others (2014) adopt a scaling law that uses the c coefficient obtained by Chen 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 Bahr and others (1997) based on physical considerations for mountain glaciers (which we note is different from the γ = 1:36 obtained by Bahr (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 Radić and others (2014) following Radić 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 Dunse (2011), Pettersson and others (2011) and Martín-Español (2013), which are 2559 km³ (Austfonna) and 442 km³ (Vestfonna), with respective relative errors of 3% and 7%, estimated by Martín-Español (2013). Thus, the ice volume of Nordaustlandet ice caps totals 3001 ± 81 km³. 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³, as our best estimate of Svalbard ice volume. In terms of sea-level equivalent (SLE), assuming an oceanic area of 3.62 × 10⁸ km² and a glacier ice density of 900 kg m⁻³, our volume estimate corresponds to a total potential contribution to sea-level rise of 17 ± 2 mm SLE.

Table 6. Estimated total glacier volume of Svalbard and its potential contribution to sea-level rise

This volume estimate is in the low range of those published in the literature (Table 1), including the regionally based estimates of Macheret and others (1984) and Hagen and others (1993), but higher than that of Grinsted (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 Macheret and others (1984) and Hagen 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 (Pfeffer 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 (Arendt 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³ (depending on the scaling relation used), which only applies to the comparison between the estimate by Marzeion 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 Huss 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 (Pfeffer 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 Grinsted, 2013). This gives a total volume of Nordaustlandet ice caps of 2444 km³ (averaging the logmse- and absdev-based results), which is 557 km³ lower than the volume of 3001 km³ obtained from the GPR data. Nearly all the difference (546 km³) corresponds to Austfonna. Using the scaling-based estimate for all of Svalbard gives a total volume of 6143 km³, making our result closer to that of Grinsted (2013). It would be even closer (6003 km³) 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³, 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 Grinsted (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³ 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 Farinotti 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.