Svalbard glacier elevation changes and contribution to sea level rise

[1] We compare satellite altimetry from the Ice, Cloud, and Land Elevation Satellite (ICESat, 2003–2007) to older topographic maps and digital elevation models (1965– 1990) to calculate long-term elevation changes of glaciers on the Svalbard Archipelago. Results indicate significant thinning at most glacier fronts with either slight thinning or thickening in the accumulation areas, except for glaciers that surged which show thickening in the ablation area and thinning in the accumulation areas. The most negative geodetic balances occur in the south and on glaciers that have surged, while the least negative balances occur in the northeast and on glaciers in the quiescent phase of a surge cycle. Geodetic balances are related to latitude and to the dynamical behavior of the glacier. The average volume change rate over the past 40 years for Svalbard, excluding Austfonna and Kvitøya is estimated to be 9.71 ± 0.55 km yr 1 or 0.36 ± 0.02 m yr 1 w. equivalent, for an annual contribution to global sea level rise of 0.026 mm yr 1 sea level equivalent.


Introduction
[2] The most recent IPCC assessment estimates that glaciers and ice caps outside Greenland and Antarctica contain between 15 and 37 cm of sea level equivalent (SLE) [Lemke et al., 2007].Even though this is small compared to the >60 m SLE of Antarctica and Greenland, it is the smaller glaciers and ice caps that are expected to be the greatest contributors to near-future sea level rise [Meier et al., 2007].Recent studies estimate that their contribution to sea level rise has been accelerating from about 0.35 -0.40 mm yr À1 SLE for the period 1960 -1990 to about 0.8-1.0mm yr À1 SLE for 2001-2004[Kaser et al., 2006]], about one third of the total observed global sea level rise.It is therefore important to quantify glacier volume changes for the various glaciated regions in the world, both to estimate glacial sea level contribution and to link such contributions to regional climatic changes.In this paper we estimate the contribution of Svalbard glaciers to sea level rise.
[3] Various methods exist to estimate regional volume changes of ice masses around the world.Traditional glacier mass balance measurements are typically extrapolated to estimate regional mass balances [Dowdeswell et al., 1997;Dyurgerov and Meier, 1997;Haeberli et al., 2007;Hagen et al., 2003aHagen et al., , 2003b]].Using mass balance data, the contribution of Svalbard glaciers to sea level rise has been estimated previously to be 0.01 mm yr À1 SLE [Hagen et al., 2003b], 0.038 mm yr À1 SLE [Hagen et al., 2003a], and 0.056 mm yr À1 SLE [Dowdeswell et al., 1997].The differences in these estimates arise from the procedures used to extrapolate traditional mass balance measurements over unmeasured areas.Hagen et al. [2003a] derive a single relation between mass balance and elevation, which is then integrated over the entire archipelago, whereas Hagen et al. [2003b] integrate 13 regional mass balance curves over the archipelago.Dowdeswell et al. [1997] use an averaged net mass balance estimated from three glaciers to integrate over the glacier area.The large variation in previous SLE estimates of Svalbard exemplifies the uncertainty in extrapolations of traditional mass balance measurements in a region where climatic spatial variability is significant.
[4] Remote sensing provides an independent approach for mass balance estimation through measurements of elevation changes using for example photogrammetry [Cox and March, 2004;Krimmel, 1999] or altimetry [Arendt et al., 2002;Howat et al., 2008;Zwally et al., 2005].Airborne laser altimetry conducted over Svalbard in 1996 and 2002 along $1000 km of profiles was too spatially limited to allow integration of the elevation changes into volume changes; however, the data suggest that eastern parts of Svalbard may be closer to mass balance equilibrium than the western and southern regions [Bamber et al., 2004;Bamber et al., 2005].Long-term volume changes estimated from maps made by a variety of methods over smaller glaciers and ice fields indicate increases in the rate of loss within the last 15 years [Kohler et al., 2007;Ka ¨a ¨b, 2008;Nuth, 2007].
[5] Satellite measurements can provide accurate estimates of recent volume and mass changes.In this paper we use the NASA Geoscience Laser Altimetry System (GLAS) instruments aboard the Ice, Cloud, and Land Elevation Satellite (ICESat) [Schutz et al., 2005].The period of ICESat observations (2003 -2007) is relatively short, and it is not always possible to distinguish snowfall and mass balance variability from true climatic signals using the satellite data alone.Longer-term comparisons, important for determining present-day anomalies, must rely on comparing the modern satellite data to older topographic data.
[6] Problems encountered with the GLAS lasers forced a greatly curtailed measurement program, both spatially and temporally.There is nevertheless sufficient data over the entire Svalbard archipelago in a 4 year period (2003 -2007) to allow comparison of ICESat elevations with older photogrammetric maps and digital elevation models (DEMs) from 1965 to 1990.This comparison is used to generate a long-term estimate of glacier volume change of various regions and subregions for the entire archipelago except Austfonna and Kvitøya.

Geographic Setting
[7] Svalbard is an Arctic archipelago ($78°N to $15°E) situated north of Norway between Greenland and Novaya Zemlya.The islands lie between the Fram Strait and the Barents Sea, which are at the outer reaches of the North Atlantic warm water current [Loeng, 1991].Therefore, Svalbard experiences a relatively warm and variable climate as compared to other regions at the same latitude.To the north lies the Arctic Ocean where winter sea ice cover limits moisture supply.To the south is a region where cyclones gain strength as storms move northward [Tsukernik et al., 2007].These geographical and meteorological conditions make the climate of Svalbard not only extremely variable (spatially and temporally), but also sensitive to deviations in both the heat transport from the south and to the sea ice conditions to the north [Isaksson et al., 2005].
[9] In this study, we divide Svalbard into five major regions; South Spitsbergen (SS), Northeast Spitsbergen (NE), Northwest Spitsbergen (NW), the Eastern Islands (EI), and Vestfonna (VF).This division derives partly from natural climatic conditions and partly from the temporal distribution of the available DEMs.In addition, subregions are defined within each region which are based upon drainage basins and the availability of spatially representative ICESat profiles.Throughout this study, two-letter abbreviations are used within the text to identify the five large regions.Three-letter codes are abbreviations for the defined subregions in the maps and tables though full names are used in the text.

Data
[10] Digitized 1:100,000 scale topographic maps made from vertical aerial photographs taken between 1965 and 1990 at scales between 1:15,000 and 1:50,000 (the Norwegian Polar Institute (NPI) S100 Series Topographic maps of Svalbard) form the base data that are compared to ICESat.In NE, NW, EI, and VF contour maps were constructed by NPI on analog stereoplotters using 1965, 1966, 1971, and 1990 imagery, respectively.The DEM for SS was constructed by NPI using the digital photogrammetry software package SOCET SET 1 , from 1:50,000 scale vertical photographs taken in 1990.The grid spacing is 20 m.Table 1 lists the regions and time intervals from which elevation changes are calculated.Austfonna and Kvitøya ice caps are not included in this analysis because the available topographic maps are of too low accuracy due to limited ground control available and due to the large low-contrast zones in the firn and snow areas which render photogrammetric elevationparallax measurement very inaccurate or impossible.The 2002-2008 volume change of Austfonna has been estimated in a separate study [Moholdt et al., 2009].
[11] ICESat contains a laser altimeter system (GLAS) that has been acquiring data since 2003.GLAS retrieves average surface elevations within $70 m diameter footprints every $170 m along track.The single shot elevation accuracy is reported to be $15 cm over flat terrain [Zwally et al., 2002], although accuracies better than 5 cm have been achieved under optimal conditions [Fricker et al., 2005].However, some data are lost to cloud cover, and ICESat performance degrades over sloping terrain and under conditions of pronounced atmospheric forward scattering and detector saturation.When the GLAS laser is transmitting pulses with high energies (i.e., during the early stages of an instrument's life) toward flat ice terrain, higher than normal echo-return energies cause detector saturation (i.e., pulse distortion) [Abshire et al., 2005].A saturation range correction [Fricker et al., 2005] available since ICESat Release 28 has been added to the elevations to account for the delay of the pulse center in saturated returns.In this study, we use the GLA06 product between 2003 and 2007 from ICESat data release 428 available from the National Snow and Ice Data Center (NSIDC [Zwally et al., 2008]).

Intersections Between ICESat Points and DEMs
[12] The satellite data, digital topographic maps, and photogrammetric DEMs are first established in the same horizontal and vertical datum and projection.The early NPI maps are referenced to the European Datum 1950 (ED50) while the 1990 DEMs are referenced to the World Geodetic System of 1984 (WGS84).A seven-parameter transformation between the UTM projections in both datums is used to convert the early maps from ED50 to WGS84.Elevation conversions are not required since the topographic maps, and photogrammetric DEMs are referenced to mean sea level using NPI mean sea level reference markers positioned around Svalbard.ICESat elevations, on the other hand, are first converted from TOPEX/Poseidon to WGS84 ellipsoid heights, and then converted to orthometric heights by subtracting the EGM96 model of geoid heights in the mean tide system.A horizontal transformation between WGS84 and TOPEX/Poseidon was not necessary since the displacements are only a few centimeters.
[13] There are various ways to produce elevation changes between ICESat profiles and contour maps.Ka ¨a ¨b [2008], for example, averages eight different methods of comparison between contours, a stereo satellite-derived DEM and ICESat profiles to estimate volume changes on Edgeøya, eastern Svalbard.For region SS, a photogrammetric raster DEM is the original data product such that elevation changes are simply calculated as differences between the ICESat point elevations and the bilinear interpolation of the underlying DEM at the locations of the ICESat footprint center.
[14] For NW, NE, EI, and VF, 50 m contours were digitized by NPI from the original map foils.Three methods for calculating the vertical differences between the contours and the ICESat-derived elevations are implemented: [15] 1. Use only ICESat points where the waveform footprint directly overlays a contour.The elevation difference is then, without any interpolation, directly calculated between the ICESat point elevation and the contour elevation included in the footprint.This method results in a small number of differences but avoids interpolation artifacts.
[16] 2. Interpolate the intersection between two successive ICESat points and a contour between them.This method results in a larger number of differences but assumes a linear slope between two successive ICESat points, i.e., over 170 m across the contours.
[17] 3. Interpolate a DEM (50 m grid spacing) from the contours using an iterative finite difference interpolation technique [Hutchinson, 1989], and subtract the DEM from the ICESat points as described above for the SS region.This method results in the largest number of differences but involves DEM-interpolation artifacts, in particular where contour lines are scarce due to low slopes.
[18] As a first measure to assess the characteristics and uncertainties of these three methods, elevation differences on nonglacier terrain, assumed to be stable, are analyzed for each region (Figure 1).The sample size of the three methods within each region is normalized to ensure proper inter-method comparison.The regional sizes are 4250 (NW), 2671 (NE), 1261 (EI), 954 (VF), and 5904 (SS) points.In all regions, method 1 results in a larger root-mean-square error (RMSE) than method 2 because elevation errors on steep slopes increase with distance between the ICESat center point and the contour.At 35 m distance, the radius of an ICESat footprint, the potential elevation error is 5 m for a 10°slope and up to 30 m for a 40°slope.The RMSE difference between methods 1 and 2 is greatest in NW, NE, and SS, where topography is dominated by jagged mountains with steep flanks rather than the plateau-type terrain characteristic as found for EI and VF.The RMSE for method 3 is significantly greater than for methods 1 and 2 in EI and VF.DEM interpolators are less accurate on terrain with large roughness (e.g., cliffs and plateau edges) and where distances between contours are large (relatively flat terrain, e.g., plateaus and strand flats).Both these topographic characteristics are predominant in EI and VF.The RMSE from method 3 is similar to that of method 2 in NW and NE since the DEM interpolation is as accurate as a linear interpolation between ICESat points in the more alpine mountainous terrain, with its dense contours and evenly steep slopes.The RMSE for method 3 in region SS is exceptionally small because the underlying raster DEM was directly measured using digital photogrammetry and did not have to be interpolated from contours.
[19] Method 2 is considered the most precise of the three methods for comparing ICESat to contours, especially over large flat surfaces.However, the distribution of ICESat-Figure 1. Box plots of the elevation differences between topographic DEMs and ICESat elevations on nonglacier terrain for the three methods outlined in section 4.1 for each region.The central point is the median, the box edges are the 25th and 75th quantile of the data, and the edges of the whiskers contain 99.3% of the data.The numbers at the bottoms of the box plots are the RMSE values of the elevation differences.contour intersections is limited in areas where profiles are transverse to the glacier centerline (for example, Figure 2a).Also, DEM interpolation of slope surfaces with limited topographic roughness, such as glaciers, introduces smaller errors than, for example, mountainous terrain outside these glaciers.Therefore, both methods 2 and 3 are implemented at the regional scale, whereas only method 3 is implemented at the subregional scale because the enhanced spatial distribution and number of elevation differences provide enough information to estimate volume change.
[20] Vertical uncertainty of the DEMs is estimated by the RMSE between ICESat elevations and the DEMs over stable terrain, assuming that the ICESat data have no error (method 3).The RMSE is largest in NW and smallest in SS (Figure 1).Accuracy of contours in the firn area may be poorer due to low optical contrast and fewer control points for the photogrammetric compilation.However, glaciers and ice caps have smoother slopes than nonglacier areas, reducing vertical errors caused by horizontal distortions and DEM interpolation.

Estimation of Elevation Change and Volume Change
[21] ICESat surface point elevations (h 1 ) are differenced to the underlying DEM pixels (h 0 ) using bilinear interpolation producing elevation changes: dh = h 1 À h 0 .Because the ICESat points are acquired in multiple years (2003 -2007), the elevation change points are divided by their respective time interval to produce point elevation change rates (dh/dt).Some outliers are present due to noisy ICESat points from atmospheric contamination, erroneous DEM elevations, or from extreme changes due to glacier surges.Outliers are removed regionally with an iterative 3s filter within 50 m elevation bins until the improvement of the resulting standard deviation (s) is less than 2% [Brenner et al., 2007].Because of variable cloud cover, some repeat track profiles are measured more often than others, biasing the spatial distribution of points toward those tracks that contain the most cloud-free profiles.Therefore, neighboring elevation differences within a 500 m radius are averaged to create one point for every kilometer along track.The original population of 92,811 elevation change points is reduced to 5631 points through this block smoothing.
[22] To regionalize thickness changes, relations describing the variation of dh/dt with elevation are created by fitting higher-order polynomial curves: dh z /dt = f(h 0 ).Higher-order polynomial fitting is less influenced by noise and outliers than averaging per elevation bin, while preserving the general trend of the elevation changes, especially at lower altitudes where thickness changes approach zero due to glacier retreat and debris-covered tongues [Arendt et al., 2006;Ka ¨a ¨b, 2008].Moreover, continuous curves allow us to estimate average thickness changes also at elevation intervals where little or no data are available due to the spatial distribution of ICESat profiles.The order of the polynomials is generally increased until the RMSE converges but also requires some subjective judgment as lowerorder fits can experience relatively low RMSE while still producing runaway tails at the edges of the data.At the regional scale, sixth-order polynomials were used while second-to sixth-order were used at the subregional scale.Glacier hypsometries (area-altitude distribution) for each where dh z /dt is the elevation change curve and A is the area at each elevation bin, Z.The dh z /dt curves and hypsometries for the five regions are shown in Figure 3.
[23] At the regional scale, a sufficient spatial distribution of dh/dt points allows robust estimates of dV ice /dt.We combine these five regional estimates to calculate the total contribution of Svalbard glaciers to sea level rise.At the subregional scale, it can be more difficult to obtain a spatially representative dh/dt distribution.Approximately 33% of the glaciated area in this study did not have a suitable spatial distribution to estimate subregional volume changes.Some glaciers have only a few ICESat profiles resulting in large data gaps at some elevation bins.In cases where these data gaps are greater than 3 -4 elevation bins, a straight line is used to interpolate dh/dt between higher and lower ICESat profiles.In cases where the data distribution was still too sparse, adjacent glaciers are combined to produce better dh/dt relationships.As an alternative to using polynomial elevation change curves to generate volume change rates, we could also have used the mean [Nuth et al., 2007] or median [Abdalati et al., 2004] dh/dt for each elevation bin.The differences in estimated volume changes between using mean, median, or polynomial fits are typically 4% -7%.
[24] We assume that all volume changes are of glacier ice [Bader, 1954] and multiply dV ice /dt by 0.917 (the density of ice) to obtain water equivalent volume change rates (dV water /dt).This assumption is valid in the ablation areas, but is more uncertain in the accumulation areas, where firn thickness or density may increase or decrease.Dividing dV water /dt by the average of new and old glacier areas [Arendt et al., 2002] provides geodetic mass balance rates ( _ b in m yr À1 ).Updated glacier areas are not yet available for the ICESat epoch, and thus we divide by the older glacier area (1966, 1971, or 1990 depending on the region) which slightly underestimates geodetic mass balances due to glacier retreat.Thickness changes (dh/dt) are given in meters of ice equivalent (m ice), while volume changes and geodetic mass balances are provided in meters of water equivalent (m w. equivalent).

Errors
[25] ICESat elevations are accurate to better than 1 m [Brenner et al., 2007].Our analysis of 237 crossover points within individual ICESat observation periods (<30 days) over Svalbard glaciers yielded a standard deviation of the elevation differences of 0.6 m, for slopes <15°.Therefore, the greatest sources of error within our estimates derive from the photogrammetrically derived topographic maps and DEMs.Errors in these products typically result from low radiometric contrast in the images, and lack of availability and quality of ground control points.To estimate single point accuracies of the DEMs relative to ICESat, we use the population of elevation differences over nonglaciated surfaces (see section 4.1) with slopes similar to those of glaciers (<15°).The population sizes for the stable terrain data sets range from 6500 points in NW to 11,826 points in SS, all distributed along the ICESat tracks.The regional ICESat-DEM differences approximate Gaussian distributions with mean differences ranging from À0.4 to 2.7 m (Table 1).We attribute these biases to ground uplift (0.10-0.24 m within our measurement periods [Sato et al., 2006]), to snow cover in the ICESat observation periods (maximum 1 m), and to deviations between the EGM96 geoid model and the mean sea level references used in NPI maps (maximum 1 m).Individual vertical biases are removed from their respective regions.
[26] Individual point elevation change uncertainties (E PT ) are estimated by the root sum squares (RSS) of the uncertainties of each data set (Table 1).Image contrast in glacier firn areas is typically low, leading to a problem known as ''floating contours'' [Arendt et al., 2002].This effect has previously been accounted for by assigning accumulation area contours a vertical uncertainty of two to three times that of an ablation area contour [Adalgeirsdottir et al., 1998;Arendt et al., 2006;Nolan et al., 2005].We use a stepwise assignment of accuracies to the different surface types by assuming that the lowermost one-third of the elevation bins for each region and subregion corresponds to the ablation zone, the middle one-third to the zone around the equilibrium line altitude (ELA), and the upper one-third to the accumulation zone.Individual point elevation change uncertainties are then estimated by where z is the respective elevation bin, s DEM is the standard deviation of terrain differences on slopes less than 15°( Table 1), s ICESat is conservatively assigned to 1 m, and c is 1, 2, or 3 for the ablation, ELA, and accumulation zones, respectively.The middle one-third of elevations for each region (Figure 3) corresponds to a rough ELA map of Svalbard [Hagen et al., 2003b].The simplicity of this parameterization for c does not warrant a precise ELA location because anything above the lowest one-third elevation bins receives an uncertainty at least double the estimated elevation errors in the ablation area.Moreover, the resulting errors for each zone in Table 1 are provided as average annual rates to emphasize the reduction of error by having a longer time span between measurements.
[27] An additional error source arises from extrapolation (E EXT ) of a limited number of dh/dt points to the entire glaciated surface.E EXT represents the uncertainty about the mean elevation change rate, estimated by the spatial variability of thickness changes rates [Arendt et al., 2006;Thomas et al., 2008].We use the standard deviation of glacial dh/dt within each 50 m elevation bin as an approximation for the extrapolation error.At the subregional scale, elevation bins that have too few measurements (less than 5 dh/dt points) are set to twice the regional mean E EXT within corresponding elevation bins.
[28] Errors in volume changes and geodetic mass balances are estimated as the combination of the two error components in each elevation bin; (1) the point elevation error (E PTz ), and (2) extrapolation error (E EXT ).Elevation changes are averaged by elevation, thus errors are reduced by the square root of the number of independent measurements within each elevation bin.Both E PTz and E EXT are random so that the combined errors are summed by RSS to produce the total elevation change error at each elevation bin (z): When full spatial coverage is available, N may be represented by the number of pixels or measurements, assuming there is no spatial autocorrelation [Etzelmu ¨ller, 2000].Conservatively, we account for spatial autocorrelation and the varying distribution of ICESat profiles over subregions and regions by defining N as the number of independent ICESat profiles within each elevation bin, rather than the number of actual data points.Thus, profiles containing more than one point within an elevation bin are, for error assessment purposes, reduced to one measurement.
In our study, the total number of ICESat footprints on glaciers (92,811) is reduced by smoothing to 5631 points (see section 4.2) whereas N for the entire study area becomes 2482.Figure 4 shows an example of the elevation dependency of the error types.The standard point errors are largest at higher elevations where poor radiometric contrast makes photogrammetry difficult while the extrapolation errors are largest at the lowest elevations where spatial variability of elevation changes is greatest due to glacier retreat and differential ablation (clean ice versus dirty ice).
[29] Volume change errors (E VOL ) are then estimated by the RSS of the elevation errors (E Z ) multiplied by the area (A Z ) assuming that the errors are independent between the elevation bins (Z): Uncertainties in the glacier outlines of the DEMs are an additional source of error.We expect this error to be small at the spatial scale of our volume change estimates, and therefore do not account for it.Updated glacier outlines are not available for the ICESat epoch, so while total volume change estimates are correct, geodetic mass balances are slightly underestimated due to the prevalent glacier retreat on Svalbard during the study period [Hagen et al., 2003b].
Errors from seasonal differences between the end of summer DEMs and the multiseasonal ICESat acquisitions (in February/March, May/June, and September/October) are not more than 1 m or $ 0.02-0.13m yr À1 over the decadal measurement period.However, after accounting for the mean ICESat-DEM bias, the seasonal acquisition of ICESat introduces errors in both directions due to snow depths in acquisitions before July/August (photographic acquisitions for the topographic maps and DEMs) and additional melt that occurs after July/August thus becoming a part of our random point error.Last, the 4 year time span of ICESat

Thickness Changes
[30] Average annual elevation change rates (dh/dt) over Svalbard are shown in Figure 5.The mean observed dh/dt for Svalbard (excluding Austfonna and Kvitøya) is À0.40 m yr À1 with 95% of the data ranging between À1.65 and +0.85 m yr À1 .Frontal thinning is observed nearly everywhere except for those glaciers that have surged in the observation period.Regionally, the most negative average annual frontal thinning rates occur on SS, NW, and EI, respectively, while NE and VF have the lowest average annual frontal thinning rates (Figure 3).At higher elevations, glaciers experience only slight vertical changes, with both thinning and thickening found.Regionally, only VF experiences an average annual thickening at higher elevations while NW, NE, EI, and SS experience thinning varying between À0.15 and À0.30 m yr À1 .Extreme dh/dt occurs where the calving fronts of marine terminating glaciers changed their position or on glaciers that have surged.For example, the minimum and maximum change rates for the entire study area (À4.92 and +8.21 m yr À1 ) are found on two glaciers that have surged recently: Perseibreen which surged in 2000 [Dowdeswell and Benham, 2003] and Fridtjovbreen, which surged in the mid-1990s [Murray et al., 2003a].
[31] Northwest Spitsbergen (NW) glaciers experience widespread frontal thinning of À1 to À2 m yr À1 .The largest frontal thinning occurs on Borebreen and in Trollheimen Land.Some glaciers in NW are thinning throughout their accumulation areas (e.g., Isachsenfonna and Holtedahlfonna) while others experience significant increases at higher elevations (e.g., Kongsvegen, Borebreen, Holmstro ¨mbreen, Morabreen, and Orsabreen).On Kongsvegen, the thickening of $0.5 m yr À1 at upper elevations and thinning of À1 m yr À1 at lower elevations correspond well to previous measurements between 1991 and 2001 [Hagen et al., 2005] but is less than that measured between 1966 and 1996 [Melvold and Hagen, 1998], +1 and À2.5 m yr À1 at upper and lower elevations, respectively.The surge on Abrahamsenbreen in 1978 [Hagen et al., 1993] resulted in +80 m vertical increases at the tongue and À40 to À80 m decrease in the reservoir, while the surge of Osbornebreen in 1987 [Rolstad et al., 1997] is seen as +50 m frontal increases and À20 m decreases at higher elevations.The larger elevation changes measured using 1966 and 1990 maps [Rolstad et al., 1997] result because the latter map was made more or less at the termination of the surge, and the glacier has been rebuilding since 1990.
[32] Frontal thinning rates in NE are more moderate than in NW, ranging between À1.5 and 0 m yr À1 , with the most negative rates occurring at the calving fronts draining Lomonosovfonna, Kongsfonna, and Negribreen.Large thinning rates (À1.5 to À0.5 m yr À1 ) also occur in the upper elevations of Tunabreen and Hinlopenbreen due to recent surges.Thickening of +0.5 to +1 m yr À1 is observed at the higher elevations of Negribreen from 1966 to $2005, slightly larger than the dh/dt measurements of +0.2 to +0.5 m yr À1 between 1996 and 2002 from airborne laser scanning [Bamber et al., 2005].The upper elevations of Kongsfonna and Balderfonna ice caps have thinned by À0.2 to À0.3 m yr À1 .A ˚sgardfonna ice cap is generally thinning (À0.1 to À0.3 m yr À1 ) at higher elevations although slight thickening is observed toward the northeast.Similar patterns and magnitudes of elevation change rates were observed at high elevations (±0.10 m yr À1 ) between 1996 and 2002 [Bamber et al., 2005].Ursafonna has thinned across the top of the ice cap (À0.1 to À0.3 m yr À1 ) although +60 to +80 m frontal increases occur at the confluence between Chydeniusbreen and Polarisbreen.No surges have previously been recorded for these glaciers.On Oslobreen, the southern outlet glacier of Ursafonna, mid-elevation thickening of +0.3 to +0.5 m yr À1 is apparent both between 1966 and 2005 (this study) and from 1996 to 2002 [Bamber et al., 2005].
[33] Frontal thinning of glaciers on the EI is most similar to that of NW, ranging between À0.3 and À2.0 m yr À1 .The largest frontal thinning occurs on Edgeøyjøkulen and Digerfonna.Elevation changes of Digerfonna are similar to those reported by Ka ¨a ¨b [2008] who use the same data as in this study as well as a DEM from ASTER satellite stereo imagery.Slight thickening is observed at the higher elevations of Barentsøyjøkulen and Edgeøyjøkulen (+0.2 m yr À1 ), where many of the outlets are suggested to be surge type [Dowdeswell and Bamber, 1995].Storskavelen, a small ice cap northwest of Edgeøyjøkulen, experiences moderate thinning (0 -1 m yr À1 ) across the entire surface.
[34] Vestfonna (VF) contains the largest dh/dt variation out of all the regions, and is the only region in which significant thickening is observed.On the south side, Aldousbreen, Frazerbreen, and Idunbreen experience frontal thinning (up to À1 m yr À1 ) and upper elevation thickening (up to +1 m yr À1 ).Gimlebreen in the southwest has thinned over the entire surface.Bodleybreen surged in the late 1970s [Dowdeswell and Collin, 1990].Since 1990, the upper glacier has thinned dramatically (À1 to À2 m yr À1 ) implying another surge may have occurred or is occurring.Franklinbreen has thickened, greatest at lower elevations, consistent with a post-1990 advance reported by Sneed [2007].To the north, both Maudbreen and Sabinebreen experience slight frontal thinning and high-elevation thickening.Rijpbreen has advanced since 1990, with midelevation thinning and high-elevation thickening.In general, the greatest thickening of the Vestfonna ice cap occurs along the northern ridge.Note however that point elevation errors on VF are estimated to be as much as ±0.6 m yr À1 at lower elevations (see section 4.3), such that most of the dh/dt values are not statistically significant.
[35] In SS, frontal thinning ranges from À0.3 to À3.0 m yr À1 .High-elevation dh/dt values range from À1.0 to +0.5 m yr À1 ; however, most of SS is thinning at rates of À0.1 to À0.3 m yr À1 .On Sørkapp, all dh/dt are negative, with frontal thinning rates up to À2.5 m yr À1 .In Wedel Jarlsberg Land, all glaciers have thinned.Thinning rates on the middle elevations of Rechecherbreen (1990Rechecherbreen ( -2005) ) from À0.5 to À1.0 m yr À1 are similar to those measured between 1996 and 2002 [Bamber et al., 2005].Zawadskibreen, Polakkbreen, and Vestre Torrellbreen experienced highelevation thinning and mid-elevation thickening due to surges in an early stage [Sund et al., 2009].In Heerland, many glaciers have thinned, with the exception of glaciers that surged recently (e.g., Skobreen, Perseibreen, and Ingerbreen), and glaciers that are potentially in a quiescent phase of a surge cycle (e.g., Indrebreen, Inglefieldbreen, and Edvardbreen), which experience mid-and high-elevation increases.Glaciers in South Sabine Land have thinned at all elevations, with values ranging from À2.7 m yr À1 at the front of Elfenbeinbreen to À0.1 m yr À1 at high elevations of Fimbulisen.

Volume Changes and Geodetic Mass Balances
[36] The average volume change rate for 27,000 km 2 glaciers on Svalbard is À9.71 ± 0.53 km 3 yr À1 w equivalent, the equivalent to a geodetic balance of À0.36 ± 0.02 m yr À1 (Table 2).The most negative regional geodetic mass balance occurs in SS.However, SS is estimated over a shorter and more recent time interval (1990 -2005), compared to the other regions.We note that the mass balance for the latter period is almost twice as negative as its 1936 -1990 geodetic balance estimate [Nuth et al., 2007], consistent with Kohler et al.'s. [2007] conclusion.Northeast Spitsbergen has the least negative balance of À0.25 ± 0.02 m yr À1 while VF is the only region with a positive geodetic balance (+0.05 ± 0.17 m yr À1 ) though not significantly different from zero.Austfonna, not included in this study, has experienced interior thickening [Bamber et al., 2004] concurrent with extensive marginal thinning and retreating.Dowdeswell et al. [2008] suggest a volume change rate between À2.5 and À4.5 km 3 yr À1 , although recent altimetric measurements for the period 2002 -2008 indicate total losses on the order of À1.3 km 3 yr À1 [Moholdt et al., 2009].
[37] Table 2 provides a list of the regions and subregions with their associated volume changes, geodetic mass balances, and error estimates.The spatial variability of the subregional geodetic balances can be seen in Figure 6.As with the regional estimates, the most negative geodetic balances occur in the south and along the western and eastern coasts.Moderately negative geodetic balances occur toward NE while the subregions of Vestfonna show the most positive balances though the largest errors.Our estimate for Digerfonna in EI of À0.49 m yr À1 is similar to Ka ¨a ¨b [2008], who estimated À0.5 m yr À1 by comparing the same NPI map data to a 2002 DEM generated from ASTER stereo imagery.In NW, the most negative balances occur on the surged glaciers of Abrahamsenbreen and Osbornebreen, and in subregions such as Trollheimen Land and Albert I Land that are thinning at most elevation bins.The least negative balances occur on those glaciers suspected to be in a quiescent phase of a surge cycle that experience thickening such as Kongsvegen, Borebreen, and Holmstro ¨mbreen.Similarly in NE, the most negative balances occur on the surged glaciers Hinlopenbreen and Tunabreen, and on Lomonosovfonna which is drained by two large tidewater glaciers.The least negative geodetic mass balance occurs on Negribreen which shows significant high-elevation thickening.In SS and EI, the most negative balances occur in the south while VF is the only region that shows a mix of positive and negative balances.
[38] Hagen et al. [2003aHagen et al. [ , 2003b] ] provide an overview of the in situ mass balance measurements available around Svalbard.The time periods for such measurements vary; however, they are the only available measurements from which to compare.In NW, Midre Love ´nbreen and Austre Brøggerbreen are among the longest arctic mass balance measurement series [Hagen and Liestøl, 1990].Their mean annual net balances of À0.38 and À0.48 m yr À1 are similar to our subregional estimate for Brøgger-halvøya and Prins Karls Forland of À0.43 m yr À1 .There is a discrepancy on Kongsvegen, however, where our estimate of À0.23 m yr À1 is significantly more negative than the published in situ mass balance estimate of 0.00 to +0.04 m yr À1 [Hagen et al., 2003a[Hagen et al., , 2003b]], or including the most recent years, À0.06 m yr À1 .Kongsvegen mass balance measurements start from 1987, representing about 20 years of data; this is only about half of our measurement period, roughly 40 years.In SS, mean net balances on Hansbreen (À0.52 m yr À1 ) and Finsterwalderbreen (À0.51 m yr À1 ) are less negative than our subregional estimate for Wedel Jarlsberg Land (À0.65 m yr À1 ), which probably relates to spatial and temporal differences between the measurements.However, they correspond well to the regional SS estimate of À0.55 m yr À1 , which also includes glaciers such as Longyearbreen, Vøringbreen, Austre, and Vestre Grønnfjordbreane.In situ measurements on the latter glaciers show balances of À0.55, À0.64 À0.46, and À0.63 m yr À1 , respectively [Jania and Hagen, 1996].

Elevation Change Estimation Methods
[39] Section 3.1 outlined three methods to derive elevation changes between ICESat and contour lines.On nonglacier terrain, method 2 proved to be the most accurate, especially in plateau-type terrain of EI and VF.Method 2, however, requires that the ICESat profiles cross contour lines.The distribution of such intersections on glacier tongues is limited in Spitsbergen, where glaciers tend to flow through steep valleys (see Figure 2a for an example).The opposite is true for the rest of Svalbard because the perimeter of lower-elevation contours on ice caps is largest.Method 3 introduces errors between contour lines, in places where the distance between contours is large.On the other hand, method 3 increases the number and spatial distribution of elevation change points.In using a ''hypsometric'' approach to estimate volume changes, we assume that the sampling distribution is representative for each elevation bin.The undersampling of method 2 at the lowest elevations has a larger impact on the volume loss of retreating glaciers than the uncertainty of the interpolation at higher elevations.Geodetic balances calculated by method 2 are 10% less  2).Numbers refer to glaciers mentioned in the text without individual estimates of volume changes and geodetic mass balances.Note that elevation changes in the subregion Brøgger-halvøya/Prins Karls Forland (BKF) within NW is from 1990 to 2005.negative than method 3 in NW and SS (Table 3) because glacier tongues in these regions experience the greatest thinning rates and are situated in glacier valleys with few intersections between ICESat profiles and contours.The methods differ only slightly in NE because of the smaller frontal thinning rates (Figure 3).On EI, the geodetic balances calculated from the two methods are similar, while the difference on Vestfonna is hardly significant.
[40] At the subregional scale, the distribution of elevation change points has a much larger impact on the estimates.[1993].Subregions with multiple names indicate either that the glacier has two names, one for the upper area and one for the tongue, or that glaciers were combined to form one subregion.Footnotes c -f represent a classification of surge/quiescent/ normal glaciers.Note, however, that the classification is not strict, and some glaciers may be surge-type though not inferred here to be so (the same for glaciers in a quiescent phase).Also, some subregions classified as normal glaciers may contain surge-type glaciers (i.e., Barentsøyjøkulen) and other subregions may contain a mix of surge/nonsurge glaciers (i.e., Wedel Jarlsberg Land) and are identified as such.
b For each group, the region is listed first, followed by the subregions.

F01008
NUTH ET AL.: SVALBARD GLACIER-SEA LEVEL RISE Figure 6.Estimated geodetic mass balances (w.equivalent) from subregional basins that contained a sufficient distribution of ICESat profiles.The gray areas in the background are glaciers that are not estimated individually ($33% of the total glaciated area), but are included in the total estimate for Svalbard.The exception is Austfonna which lies to the east of Vestfonna, and is not included within this study.Also shown is the number of original ICESat elevation change points (n) resulting from each of the three methods.The variation between the methods provides a reliability assessment based on different geodetic methods to compare the ICESat points to contours.
Errors in volume changes and geodetic mass balances are therefore considerably larger for the subregions than for the regions (Table 2).The sensitivity of the polynomial fitting to these limited data sets is tested on Holtedahlfonna, where airborne laser altimetry on the centerline was acquired in 2007 by the National Space Institute at the Technical University of Denmark (Figure 2).No formal error estimate has been made on the latter data, but airborne laser altimetry typically achieves accuracies on the order of no more than a few tens of centimeters.The centerline profiles are compared to the DEM via method 3 and a polynomial is fit following the same procedures as in section 4.3.Despite only five transverse ICESat tracks across the glacier, the dh/ dt curve from ICESat is similar to that produced by the centerline profile (Figure 2b).The geodetic balance calculated from the centerline laser altimetry is À0.40 m yr À1 , $20% less than our estimate of À0.49 ± 0.13 m yr À1 but within the error estimate.
[41] Overall, the regionalization of multitemporal thickness change rates from ICESat to map comparisons provides relatively robust estimates of volume changes and geodetic mass balances as measurements from extreme mass balance years will be smoothed out in the overall change rates due to the many different time spans involved.In addition, the multiseasonal acquisition of ICESat (February/March, May/June, and September/October) introduces a quasirandom error (see section 4.3) that limits the need for seasonal elevation adjustments, typically important in areas of high melt [Andreassen et al., 2002;Cox and March, 2004].The ICESat point profiles are distributed in many orientations over many glaciers rather than being concentrated along the centerlines of selected glaciers, as in other altimetric studies [Abdalati et al., 2004;Arendt et al., 2006Arendt et al., , 2002;;Bamber et al., 2005;Echelmeyer et al., 1996].Unfavorable track configurations make it more difficult to estimate geodetic mass balances from ICESat on individual glaciers, but on a regional scale the estimate benefits from containing a greater spatial distribution of points.In addition, the location of ICESat tracks relative to the glaciers is random in regional estimates, reducing the risk of systematic errors from measurement positions.

Geodetic Balance Uncertainties
[42] The geodetic mass balance on a glacier is estimated by dividing the total volume change by the average of the old and new areas [Arendt et al., 2002], and should be equivalent to the mass balance measured in situ [Elsberg et al., 2001;Krimmel, 1999].Since updated glacier masks for the ICESat epoch are not available, we underestimate geodetic balances because most of the glaciers are retreating.Glacier areas in SS retreated by $0.3% per year between 1936 and 1990 [Nuth et al., 2007].This rate is used to coarsely predict the glacier area during the ICESat epoch.Re-calculation of geodetic balances using the new predicted area results in an average difference of $5% which we take to represent the expected underestimation of our estimates.However, the final geodetic balance estimates in Table 2 do not consider area changes.
[43] The calculation of a water equivalent geodetic mass balance and sea level contribution assumes that elevation changes are the result of changes in ice thickness rather than variations in firn thickness.It is difficult to test this assumption because firn thickness varies in time and space.Long-term records (i.e., ice cores) exist only at single points at specific times.On Holtedahlfonna in NW and Lomonosovfonna in NE, the firn thickness was measured to be $20-30 m [Kameda et al., 1993;Pohjola et al., 2002b;Sjo ¨gren et al., 2007].On A ˚sgardfonna, Vestfonna, and Austfonna, the firn thickness was measured to be only 6 -10 m [Brandt et al., 2008;Dunse et al., 2009;Pinglot et al., 2001;Uchida et al., 1993;Watanabe et al., 2001].Deep firn density curves are lacking in SS and EI.The available density curves (Figure 7) show that the depths to the firn-ice transition on Holtedahlfonna and Austfonna have not changed significantly within the 15 and 8 year time intervals, respectively.
[44] Alternatively, one can apply an exponential density function (for example, r(z) = 0.9 À k(z) a ) where k and a are tuning parameters, and z is the elevation bin.The function is forced such that the lowest elevations receive a water equivalent conversion of 0.917 where higher elevations gradually receive a lower conversion factor approaching 0.55.Re-calculation of volume changes results in a 3% -7% difference in estimates.It is likely that thinning glaciers will lose some firn as the ELAs rise.This would cause an Profiles from Lomonosovfonna (LMF), 1997 [Pohjola et al., 2002b]; Holtedahlfonna (HDF), 1992 [Uchida et al., 1993], and 2005[Sjo ¨gren et al., 2007].Firn-ice transition is between 15 and 19 m. (bottom) Profiles from Vestfonna (VF), 1995 [Watanabe et al., 2001]; Austfonna, 1999 [Pinglot et al., 2001], and 2007 [Brandt et al., 2008;Dunse et al., 2009].Firn-ice transition in Nordaustlandet ranges between 7 and 10 m.Density curves on Holtedahlfonna and Austfonna are similar despite 15 and 7 year time difference, respectively, suggesting no significant changes in the firn thickness.
overestimation in the geodetic mass balances since the density of firn is lower ($550 kg m À3 ) than that of ice (917 kg m À3 ).Setting the water equivalent conversion to that of firn in the middle one-third of elevation bins results in a maximum volume change overestimation of 15% -20%.However, in reality the loss of firn around the ELA will only be a small fraction of this.These tests exemplify that the exact water equivalent conversion will have an effect on the final water equivalent volume change but only to a small degree because at least 50%-60% of the volume changes occur in the lowest one-third of the elevation bins, that is, on ice rather than firn.
[45] In Svalbard, all marine terminating glaciers are grounded [Hagen et al., 2003b] while our elevation change measurements provide thickness change rates for ice above sea level.Therefore, thinning rates of marine terminating glaciers are underestimated within the retreat areas due to ice loss below sea level.This unaccounted mass loss is important to consider when interpreting volume changes and geodetic mass balances of tidewater glaciers at a subregion scale.Dowdeswell et al. [2008] found that the ice loss from ice-marginal retreat at Austfonna is as much as 1.4 km 3 yr À1 based on ice thickness data and retreat rates from optical imagery.Hagen et al. [2003b] assumed an average retreat rate of 10 m yr À1 and an average ice thickness of 100 m along the $1000 km long calving front of Svalbard glaciers to estimate a marine retreat loss of 1 km 3 yr À1 for the entire archipelago.The underestimation in our volume change estimates should be well below this value since we exclude Austfonna and account for marine ice losses above sea level.Applying the same assumptions to the $30 km calving fronts of Vestfonna results in a marine retreat loss on the order of 0.03 km 3 yr À1 , which is within our error estimates.
[46] When we convert the total Svalbard volume change into SLEs, the marine retreat of grounded glaciers has only a minimal effect since the ice masses below sea level are already displacing seawater.In fact, since the density of ice is less than water, ice mass loss below sea level slightly decreases Svalbard's contribution to sea level.Nonetheless, we choose not to account for this effect due to the lack of information about ice thicknesses and retreat rates of tidewater glaciers.

Interpretation of Elevation Changes
[47] In general, three geometric patterns of elevation changes are observed on Svalbard.The most predominant pattern is recognized by large frontal thinning with slight thinning at higher elevations (above the ELA).The second pattern is characterized by large frontal thinning and highelevation thickening.The third pattern is seen on glaciers that surged, a frontal thickening and high-elevation thinning.Variations to these patterns may occur when surges are still in progress at the time of the second elevation data acquisition (for example, see Sund et al. [2009]).At the lowest elevations, thickness changes are the result of ice melting and changes in ice flux.At higher elevations, Svalbard glaciers are generally thinning at rates from À0.1 to À1 m yr À1 .In Svalbard, the end of the Little Ice Age (LIA) occurred around the 1920s [Nordli and Kohler, 2004] corresponding with the onset of glacier retreat and negative mass balances [Hagen et al., 1993].Thinning above the ELA may then be explained by a decrease in the thickness of the firn column which is not the case for Holtedahlfonna (Figure 7).An alternate hypothesis may be that ice submergence velocities are larger than accumulation inputs implying a delayed or prolonged response of the dynamic system to past mass balance conditions.
[48] Some glaciers experience significant elevation increases at higher elevations.Since accumulation is dependent on atmospheric circulation and orographic effects, one would expect thickening from precipitation increases to occur regionally.Also, if accumulation is increasing, elevation increases may result from a time lag between increased mass input and the densification processes (i.e., compaction) that convert firn to ice.However, this time lag is probably shorter than the time between measurements.It is plausible that some elevation increases may be due to local precipitation increases (i.e., Asga ˚rdfonna and Vestfonna), though previous ice core research in other parts of Svalbard does not show any significant increases in accumulation rates since 1950 [Pinglot et al., 1999, Pohjola et al., 2002a].Assuming that firn density profiles remain unchanged, elevation increases above the ELA probably result from a reduced downward flux of ice that is not in balance with the present climate, most likely attributed to surge-type glaciers in quiescent phase.
[49] It remains difficult to interpret geometric changes of glaciers because a change in elevation is the result of both the mass balance and the dynamical flux [Paterson, 1994].Melvold and Hagen [1998] and Hagen et al. [2005] show that geometric changes on Kongsvegen, a surge-type glacier in quiescent phase, are equivalent to the mass balance because dynamics are stagnant after the surge in 1948.Pinglot et al. [1999] measured the mean accumulation rates from numerous ice cores around Svalbard through radioactivity measurements and dating by nuclear tests in 1963 and the Chernobyl accident in 1986.Their analysis of two ice cores on Kongsvegen results in mean net accumulation rates of +0.53 to +0.62 m yr À1 from $1963 to 1991 which compares well to our dh/dt measurements of $+0.5 m yr À1 .On Holtedahlfonna, two ice cores resulted in accumulation rate estimates between +0.47 and +0.57m yr À1 [Pinglot et al., 1999] where our dh/dt at the highest elevations show decreases of À0.25 m yr À1 .This implies a submergence ice flux of $À0.75 m yr À1 , which although quite large is not unlikely considering that Kronebreen, one of the fastest glaciers in Svalbard, drains Holtedahlfonna.Both Holtedahlfonna and Kongsvegen are situated within the same region, yet dh/dt measurements show completely different signals.Caution should be used when interpreting elevation changes, especially in a climatic context, as dynamic effects can be a major factor.
[50] A geodetic mass balance is a volume change rate normalized by the hypsometry and is thus a combination of the longer-term mass balance conditions as well as the dynamical conditions which lead to ice emergence or submergence and possibly calving.In the case of Kongsvegen above, the volume change and geodetic mass balance are solely related to the surface mass balance of the glacier since the dynamical component is essentially zero.Geodetic mass balances on surging glaciers require some care in interpreting since the changes reflect the presurge mass balance history, the ice volume lost into the sea through surging, and a postsurge mass balance history.Within this data set, the geodetic balances of surged glaciers tend to be more negative than glaciers that have not surged.A good example is that of Hinlopenbreen and Negribreen in the NE region.They are adjacent basins that have surged and are building up to a new surge, respectively.The last surge of Negribreen in 1935/1936 is reported to be one of the largest known surges in Svalbard [Hagen et al., 1993].The geodetic balance of Hinlopenbreen is the most negative in the NE region, partly reflecting the removal of $2 km 3 of ice by the 1973 surge [Liestøl, 1973] although this is not enough to completely explain the enhanced long-term volume loss.Negribreen has the least negative geodetic balance within NE due to an extensive thickening in the accumulation area that almost balances the large thinning rates on the tongue.Down-glacier transport of ice through surging should certainly lead to an immediate change in the mass balance regime by increasing the effective ablation zone, and conversely, decreasing the accumulation zone.Furthermore, crevassing increases the surface area of the ablation zone significantly, potentially also increasing melt [Muskett et al., 2003].
[51] Previous work has suggested latitude, after accounting for elevation, can explain up to 59% of elevation change variation on Svalbard [Bamber et al., 2005].In our data sets, elevation can significantly explain $30%-70% of the variation of individual dh/dt points within each region and subregion.We further test whether the volume change after normalizing by area (i.e., geodetic balance) has any spatial trends.A multiple linear regression applied to the subregional geodetic balances (population size = 37) against latitude and longitude determined that only latitude was significant in explaining 32% of the variation (a = 0.01).Removing surging (n = 5) and quiescent phase glaciers (n = 7) from the data set (population size = 25) increases the explained variance to 46% (a = 0.01).Analysis of variance (ANOVA) tests show that surged glaciers (as classified in Table 2) are more negative than quiescent-phase glaciers (a = 0.02).Variations in geodetic balances can partly be explained by latitude with southern glaciers more negative than northern glaciers and partly by the dynamical situation with surged glaciers more negative than those in a quiescent phase.

Conclusions
[52] Until now, ICESat has mainly been applied to ice sheet terrain in Antarctica and Greenland.Here, ICESat laser altimetry proved to be a highly valuable data set for estimating the regional-scale glacier volume changes for smaller glaciers and ice caps at high latitudes and with mountainous topography.The precision of ICESat elevations provides good ground control on DEM generation from satellite imagery [Berthier and Toutin, 2008;Korona et al., 2009] but also for determining the uncertainty of older topographic maps.ICESat's applicability for smaller glaciers in mountainous regions is limited because the spatial distribution of tracks does not necessarily lie along the glacier centerlines, as would be the case with airborne laser altimetric profiles.However, the spatial distribution of ICESat tracks in Svalbard is sufficient over larger regions to estimate dh/dt variation with elevation, assuming dh/dt contains normal distributions within elevation bins.The annual average volume change estimates are relatively robust due to the long time span (15 -39 years) and to the large number of measurements that 4 years of ICESat tracks provide.Errors associated with such estimates are smaller at the regional scale than at the subregional scale, mainly because of the sampling distribution.Overall, ICESat has proved to be a valuable tool to measure glacier elevation and volume changes in the Svalbard archipelago.
[53] Surface elevation changes on Svalbard glaciers vary largest with elevation, before latitude.In general, glaciers are thinning at lower elevations except on glaciers which have recently surged, where thickening is observed.At higher elevations, three change patterns are present.On some glaciers slight thinning (<À0.5 m yr À1 ) occurs in the uppermost areas which may reflect a delayed dynamic adjustment to past mass balance conditions.On other glaciers, upper altitudes are thickening, which we generally attribute to build up in the quiescent phase, as the thickening tends to occur on glaciers that have surged in the past (i.e., Negribreen and Kongsvegen), and at a subregional scale rather than at a regional scale (i.e., Borebreen and Indrebreen).Last, some glaciers experience large thinning (>À1 m yr À1 ) at the upper altitudes (i.e., Hinlopenbreen, Ingerbreen, and Polakkbreen) implying surge activity between the measurements.
[54] Geodetic balances are a measure of the long-term total glacier change, and thus reflect general climatic influences as well as local dynamic effects, mainly in surge-type glaciers.Surged glaciers have a more negative geodetic balance than neighboring glaciers that have not surged.Glaciers that seem to be in a quiescent phase of a surge-cycle have less negative geodetic balances than other glaciers.Ignoring surge-type glaciers, 46% of the geodetic mass balance variation can be explained by latitude.Spatially, the most negative geodetic balances occur in SS followed by the coastal regions of NW, Edgeøya, and Barentsøya.The glaciers in NE show moderate losses while Vestfonna is observed to be close to zero mass balance due to a moderate interior thickening that balances frontal thinning.
[55] In summary, the total volume change for Svalbard glaciers (excluding Austfonna and Kvitøya ice caps) over the past 15-40 years is À9.71 ± 0.53 km 3 yr À1 or À0.36 ± 0.02 m yr À1 w. equivalent.This corresponds to a global sea level rise of about +0.026 mm yr À1 SLE, a value which lies between two previous estimates (+0.01 and +0.038 mm yr À1 SLE) of Svalbard's contribution to sea level rise over the past 40 years [Hagen et al., 2003a[Hagen et al., , 2003b]].Our estimate is about half of the SLE contribution as estimated by Dowdeswell et al. [1997], and about 85% of the SLE contribution from the estimated volume changes published by Dyurgerov and Meier [2005].Gravity observations from the Gravity Recovery and Climate Experiment (GRACE) satellite mission between 2003 and 2007 indicate mass losses of 8.8 Gtons yr À1 [Wouters et al., 2008] which is similar to our 9.71 Gton yr À1 estimate despite the different time periods of the studies.Globally, the Svalbard contribution to sea level rise is about 4% of the total contribution from smaller glaciers and ice caps, which roughly corresponds to the area ratio between Svalbard glaciers and the sum of global glaciers and ice caps [Kaser et al., 2006].The average annual volume loss from Svalbard is about twice the 1952 -2001 loss rates in the Russian Arctic [Glazovsky and Macheret, 2006;Meier et al., 2007] and about 40% of the 1995 -2000 loss rates in the Canadian Arctic [Abdalati et al., 2004].When glacier area is considered, the Svalbard geodetic mass balance is the most negative in the Arctic, twice as negative as the Canadian Arctic, and almost four times as negative as the Russian Arctic.Lower latitude glacier regions such as Alaska [Arendt et al., 2006], Iceland [Bjo ¨rnsson and Pa ´lsson, 2008], and Patagonia [Rignot et al., 2003] are losing ice at a faster rate than Svalbard.
[56] Acknowledgments.We would like to thank the Associate Editor, Gordon Hamilton, whose comments and suggestions considerably improved the content and presentation of this article.We would also like to thank Hester Jiskoot and two anonymous reviewers for their thorough and constructive comments.This research was supported by the IPY-GLACIODYN project (176076) funded by the Norwegian Research Council (NFR).We would like to acknowledge the Norwegian Polar Institute Mapping Department for providing the maps and DEMs, and NASA and NSIDC for providing the ICESat which is a remarkably accurate data set.The National Space Institute at the Technical University of Denmark provided the 2007 airborne laser altimetry profile over Holtedahlfonna.We wish to thank Elisabeth Isaksson and Veijo Pohjola for providing firn density curves.We would also like to thank Thomas Vikhamar Schuler for improving this manuscript and providing valuable insight and discussion on the validity of ''Sorge's law.''

Figure 2 .
Figure 2. (a) Holtedahlfonna (orange basin outline) and Isachsenfonna (brown basin) in NW.Location of ICESat points are shown with dh/dt values indicated by color and method used by symbol.Black lines show the 2007 centerline airborne laser altimetry profile (The National Space Institute at the Technical University of Denmark).Background is an ASTER image from 12 July 2002.(b) Holtedahlfonna dh/dt points measured from ICESat and from the centerline laser altimetry along with their corresponding polynomial relationships with elevation.

Figure 3 .
Figure3.Area/altitude distributions of the five regions (gray, scale to the left) and the chosen polynomial fit to dh/dt by elevation (black bold line, scale to the right).The lighter black lines are one standard deviation of the original data points from the median to show the spread of the data.The lighter gray line is the number of smoothed ICESat points (see section 4.2, 5631 points for the entire study area) within each elevation bin multiplied by two for scaling purposes.The number of points thus corresponds to the y axis on the left divided by two.

Figure 4 .
Figure 4.The estimated standard point (E PT), extrapolation (E EXT ), and combined errors (E Z ) here exemplified for region NE.E PT increases with elevation because higher elevations have poorer image contrast, and thus less accurate contour elevations from photogrammetry.E EXT increases toward lower elevations because the spatial variability of dh/dt is larger because of differential ablation.

Figure 5 .
Figure 5. Annual elevation change rates (dh/dt) in (a) NW, (b) VF, (c) NE, (d) SS, and (e) EI obtained by comparing ICESat profiles from 2003 to 2007 to older DEMs from 1965 to 1990.The maps are projected in WGS84-UTM33X.Twoletter abbreviations represent the five larger regions, three-letter codes are the smaller subregions (Table2).Numbers refer to glaciers mentioned in the text without individual estimates of volume changes and geodetic mass balances.Note that elevation changes in the subregion Brøgger-halvøya/Prins Karls Forland (BKF) within NW is from 1990 to 2005.

Table 1 .
Data Sources, Time of Acquisition, Bias From Stable Terrain, and Estimated Errors for Each Region a a See equation (2) for error definitions.Individual point errors are defined for the ablation area, for the area around the ELA and the firn area.

Table 2 .
Regional and Subregional Areas, Volume Changes, and Geodetic Mass Balances in Water EquivalentAlong With Error Estimates a Abbreviations correspond with those in Figures5 and 6.Glacier names correspond to those published byHagen et al. a

Table 3 .
Regional Geodetic Mass Balance Estimates in Water Equivalent as Estimated From Using the Three Methods Described in Section 4.1 a