Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry

Remote Sensing of Environment 114 (2010) 2756–2767

Contents lists available at ScienceDirect

Remote Sensing of Environment j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / r s e

Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry Geir Moholdt a,⁎, Christopher Nuth a, Jon Ove Hagen a, Jack Kohler b a b

Department of Geosciences, University of Oslo, Box 1047 Blindern, NO-0316 Oslo, Norway Norwegian Polar Institute, Polar Centre, NO-9296 Tromsø, Norway

a r t i c l e

i n f o

Article history: Received 22 March 2010 Received in revised form 18 June 2010 Accepted 21 June 2010 Keywords: ICESat Glaciers Ice caps Svalbard Elevation changes Volume changes Mass balance Sea level change Laser altimetry

a b s t r a c t We have tested three methods for estimating 2003–2008 elevation changes of Svalbard glaciers from multitemporal ICESat laser altimetry: (a) linear interpolation of crossover points between ascending and descending tracks, (b) projection of near repeat-tracks onto common locations using Digital Elevation Models (DEMs), and (c) least-squares fitting of rigid planes to segments of repeat-track data assuming a constant elevation change rate. The two repeat-track methods yield similar results and compare well to the more accurate, but sparsely sampled, crossover points. Most glacier regions in Svalbard have experienced low-elevation thinning combined with high-elevation balance or thickening during 2003–2008. The geodetic mass balance (excluding calving front retreat or advance) of Svalbard's 34,600 km2 glaciers is estimated to be −4.3 ± 1.4 Gt y−1, corresponding to an area-averaged water equivalent (w.e.) balance of −0.12 ± 0.04 m w.e. y−1. The largest ice losses have occurred in the west and south, while northeastern Spitsbergen and the Austfonna ice cap have gained mass. Winter and summer elevation changes derived from the same methods indicate that the spatial gradient in mass balance is mainly due to a larger summer season thinning in the west and the south than in the northeast. Our findings are consistent with in-situ mass balance measurements from the same period, confirming that repeat-track satellite altimetry can be a valuable tool for monitoring short term elevation changes of Arctic glaciers. © 2010 Elsevier Inc. All rights reserved.

1. Introduction Satellite radar altimetry has been used to measure elevation changes in Greenland and Antarctica since the late 1970s (e.g. Zwally et al., 1989; Wingham et al., 1998; Johannessen et al., 2005). The large footprint size of satellite altimeters has made it difficult to apply these measurements to higher relief glaciers and ice caps. However, newer, higher resolution altimeters like the CryoSat-2 radar altimeter (Wingham et al., 2006) and the ICESat laser altimeter (Zwally et al., 2002) provide elevation data sets that can be compared to maps/ DEMs (e.g. Sauber et al., 2005; Muskett et al., 2008; Nuth et al., 2010), to airborne altimetry (e.g. Thomas et al., 2005) and to each other (e.g. Smith et al., 2005). The most established technique to obtain elevation changes directly from satellite altimetry is to compare elevations at crossover points between ascending and descending satellite passes. This is a very accurate method (Brenner et al., 2007), but the spatial sampling is typically too coarse for volume change calculations apart from in Greenland and Antarctica. Repeat-track analysis provides a much denser sample of elevation change points, but sacrifices accuracy due to the imprecise repetition of satellite ground tracks. Still, ICESat near repeat data have been used to identify grounding zones of ice shelves (Fricker and Padman, 2006), to map subglacial lakes and drainage (Fricker et al., 2007; Smith et al., 2009), and to ⁎ Corresponding author. Tel.: + 47 99102900. E-mail address: [email protected] (G. Moholdt). 0034-4257/$ – see front matter © 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2010.06.008

quantify elevation change rates in Greenland and Antarctica (Howat et al., 2008; Slobbe et al., 2008; Pritchard et al., 2009). Arctic glaciers and ice caps are among the largest contributors to sea level rise (Kaser et al., 2006). In-situ mass balance measurements are sparse in these regions, implying a need for remote sensing data to better understand regional variations in mass balance. The most used techniques to obtain elevation changes in the Arctic have been to compare multi-temporal photogrammetric maps/DEMs (e.g. Nuth et al., 2007; Kääb, 2008) or repeated airborne laser profiles (Abdalati et al., 2004; Bamber et al., 2005). However, airborne campaigns are expensive, and photogrammetry is difficult in the accumulation areas of large ice caps where there are few ground control points and where the image contrast is poor. ICESat altimetry data are freely accessible (Zwally et al., 2008) and provide a dense spatial and temporal coverage of high quality elevation points in these high latitude regions. In this article, we investigate the potential of repeat-track ICESat altimetry to derive short term glacier elevation changes within a semi-alpine high latitude environment like the Svalbard archipelago in the Norwegian Arctic. Two methods of repeat-track analysis are tested, and the results are validated against crossover points and external DEMs. Area-averaged 2003–2008 elevation change rates are estimated for 7 glacier regions as well as for the entire archipelago. Additionally, ICESat's 2–3 observation campaigns per year provide the opportunity to calculate winter and summer elevation changes. The area-averaged seasonal estimates are compared and validated with surface mass balance data from the same period.

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

2. Svalbard glaciers Svalbard is an Arctic archipelago located north of Norway between Greenland and Franz Josef Land. There are six major islands (Fig. 1), of which Spitsbergen is the largest one and the only one with permanent settlements. The meteorological conditions are temporally and spatially diverse due to the location at the confluence zone between cold and dry polar air masses from the north and more warm and humid air masses from the Atlantic currents to the southwest. Rainfall and snowfall can happen at any time of the year, and the temperature fluctuations are large, especially during winter. The large year-to-year variations in seasonal temperatures and precipitation imply that climatic trends need to be strong or averaged over long time series in order to be statistically significant (Førland & Hanssen-Bauer, 2003). The total glaciated area on Svalbard is 34,560 km2, representing about 6% of the worldwide glacier cover outside of Greenland and Antarctica. Spitsbergen is the most alpine island, having small cirque glaciers as well as extensive ice fields and valley glaciers. The eastern islands, facing the Barents Sea, have less relief and are dominated by low-elevation ice caps. Most glaciers and ice caps are considered to be polythermal (e.g. Bjornsson et al., 1996), and 60% of the glaciated areas drain into tidewater glaciers (Blaszczyk et al., 2009). Glacier dynamics are typically slow with velocities b10 m y−1 (Hagen et al.,

2757

2003a), but surge activity has been observed over most of Svalbard (e.g. Lefauconnier & Hagen, 1991; Hamilton & Dowdeswell, 1996), also in recent years (e.g. Sund et al., 2009). Annual mass balance records of small glaciers in western Spitsbergen indicate a negative mass balance regime since at least the mid 1960s (Hagen et al., 2003b). Comparisons of photogrammetric maps/DEMs, dating back to 1936, show substantial decreases of glacier area and volume (Nuth et al., 2007) with enhanced thinning rates after 1990 when compared to recent airborne lidar (Bamber et al., 2005; Kohler et al., 2007) and ICESat altimetry (Nuth et al., 2010). The mass balance of northeastern Spitsbergen glaciers has been less negative than the western ones (Bamber et al., 2005; Nuth et al., 2010). The Nordaustlandet ice caps, Austfonna and Vestfonna, have been close to balance over the last two decades (Pinglot et al., 2001; Moholdt et al., 2010; Nuth et al., 2010) if the calving front retreat losses are ignored (Dowdeswell et al., 2008). We have divided Svalbard into seven glacier regions (Fig. 1): Northwestern Spitsbergen (NW), Northeastern Spitsbergen (NE), Southern Spitsbergen (SS), Barentsøya and Edgeøya (BE), Vestfonna ice cap (VF), Austfonna ice cap (AF) and Kvitøyjøkulen ice cap (KJ). The regions are mostly consistent with Nuth et al. (2010), but we have added Austfonna (AF) and Kvitøyjøkulen (KJ) ice caps, while the Southern Spitsbergen region (SS) is restricted

Fig. 1. Average 2003–2008 glacier elevation change rates (dh/dt) across the Svalbard archipelago. Most of the dh/dt rectangles represent the plane method, but clusters from the DEM method are also present where there are no planes. The 7 glacier regions (Table 1) used in the analysis are also outlined.

2758

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

to the glaciers south of the Nordenskiöld peninsula (NP) which has too little data to be included.

data from the summer campaigns, reflecting the meteorological conditions in Svalbard with more cloud cover in summer than winter.

3. Data 3.2. Glacier DEMs 3.1. ICESat laser altimetry The Geoscience Laser Altimeter System (GLAS) onboard ICESat has been operating over 15 observation campaigns of ~35 days between 2003 and 2008 (Fig. 2). GLAS derives ranges from the time delay between 1064 nm laser pulse transmissions and surface echo returns (Zwally et al., 2002). The ground footprints are spaced at 172 m alongtrack and have a varying elliptical shape with average dimensions of 52 × 95 m for Laser 1 and 2 (until summer 2004) and 47 × 61 m for Laser 3 (since fall 2004) (Abshire et al., 2005). GLAS was designed to achieve a single shot elevation accuracy of 0.15 m over gently sloping terrain (Zwally et al., 2002), but accuracies better than 0.05 m have been demonstrated under optimal conditions (Fricker et al., 2005). However, the performance degrades over sloping terrain (Brenner et al., 2007) and under conditions favourable to atmospheric forward scattering and detector saturation (Fricker et al., 2005). We used the GLA06 altimetry product release 28 (Zwally et al., 2008) which is based on the ice sheet waveform parameterization. A saturation range correction available in release 28 was added to the elevations to account for the delay of the pulse center in saturated returns. Each observation campaign since fall 2003 contains elevation data from a 33-day sub-cycle of the nominal 91-day repeat orbits. The satellite orbits are maintained to follow a set of ideal reference ground tracks (Schutz et al., 2005). Typical reference tracks at Svalbard contain data from 5 to 10 profiles within a ground swath of a few hundred meters. Signal absorption in optically thick clouds causes some profiles to be incomplete, while other profiles are entirely lacking. The total number of measurements varies considerably between the three annual observation campaigns (two since 2006) in Feb./Mar., May/Jun. and Oct./Nov. (Fig. 2). There are most data from the winter campaigns and least

Digital elevation models (DEMs) were used to project ICESat repeat-tracks onto common locations and to extrapolate elevation changes to unmeasured areas. The data source for DEMs and glacier outlines in the regions BE, VF and KV were 1:100000 topographic maps constructed from vertical aerial photos by the Norwegian Polar Institute. Continuous 50 × 50 m DEMs were generated using an iterative finite-difference interpolation technique (Hutchinson, 1989) on digitized data points along 50 m contours (Nuth et al., 2010). This method provides smooth DEMs without the terraced effect that can result from other interpolation techniques using contour data (Wise, 2000). The glacier DEMs were based on imagery from 1971 (BE), 1977 (KV) and 1990 (VF). Vertical root-mean-square (RMS) errors of ~10 m were estimated from comparisons with ICESat points over non-glacier terrain in slopes b15˚. In the Spitsbergen regions (NW, NE and SS) we used new 2007/ 2008 DEMs of 40 m pixel resolution from the IPY SPIRIT project (Korona et al., 2009). SPIRIT DEMs are generated from high resolution along-track SPOT 5 HRS stereoscopic images using an automatic processing scheme that rely on the orbital positioning of the satellite rather than ground truth. The optical contrast of the Spitsbergen images was mostly good, limiting the need for interpolation in uncorrelated areas where the image matching failed (Berthier & Toutin, 2008). We co-registered the DEMs to ICESat by minimizing the cosinusoidal dependency between aspect and vertical deviation over land (Kääb, 2005). The DEM in SS was made entirely from a cloud free September 1 2008 image pair, while the DEMs in NW and NE were combined from four individual DEMs generated from image pairs acquired between June and September in 2007 and in 2008. We mosaiced the four DEMs in such a way that cloudy areas with low image correlation were not used in the final DEM. The DEM accuracy was checked against the February 2008 ICESat campaign, yielding RMS elevation differences within ±5 m both on ground and on glacier ice. New glacier outlines were manually digitized from the same 2007/ 2008 imagery. About 1/4 of the NE region had to be filled in with data from ~ 1966 topographic maps since no SPIRIT products were available there. Photogrammetric mapping at Austfonna ice cap (AF) is difficult due to the extensive featureless landscape. We used a new 50 × 50 m DEM constructed from differential SAR interferometry using ICESat altimetry as ground reference. The ICESat points were used to refine the interferometric baseline, i.e. to reconstruct the geometry of the SAR acquisitions. Hence, the local slopes of the DEM around the ICESat tracks should not be affected. The RMS error of the DEM was 13 m as compared to the same ICESat profiles and 7 m as compared to independent GNSS surface profiles in the interior of Austfonna. New glacier outlines for Austfonna and the surrounding smaller ice caps were digitized from a SPOT 2008 scene for the northern and western parts and a Landsat 2001 scene for the southeastern coast. 3.3. In-situ mass balance data

Fig. 2. ICESat observation campaigns and the temporal distribution of data over glacier terrain in Svalbard. Grey bars show the total number of ICESat footprints in each campaign, while the colored bars show the number of footprints used in the dh/dt calculations for the crossover method (green), the DEM method (red) and the plane method (blue). The crossover data set is upscaled by a factor of 10 to make the bars visible. Each observation campaign spans ~35 days within the months Feb./Mar. (winter), May/Jun. (summer) and Oct./Nov. (fall). The early 2003 data were not used since the 8-day repeat orbits have not been repeated afterwards.

We used 2003–2008 in-situ mass balance data from Kongsvegen in Northwestern Spitsbergen, Hansbreen in Southern Spitsbergen, and Etonbreen at Austfonna to compare with calculated winter and summer elevation changes of Svalbard glaciers (Fig. 3). Winter mass balances are obtained in late April/early May from snow depth soundings, snow pit density measurements and stake height measurements relative to the previous summer surface, while summer mass balances derive from stake height changes and firn densities at the end of the summer, typically in late August or early

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

2759

ablation period (2–4 months). The ICESat summer campaigns in May/ June 2004–2006 were not considered in the seasonal analysis since the amount of data was lower than for the other campaigns (Fig. 2). 4.1. Elevation changes at crossover points

Fig. 3. Cumulative area-averaged elevation changes (dh) for the 2003–2008 winter seasons (October–March) and summer seasons (March–October) for all Svalbard glaciers. For comparison, cumulative mass balance curves are included for Kongsvegen (Northwestern Spitsbergen), Hansbreen (Southern Spitsbergen) and Etonbreen (Austfonna). Note that the seasons of elevation change and mass balance are spanning different time intervals (Oct.↔Mar. vs. Sept.↔May), and that the area-averaged elevation changes (m) differ from the water equivalent mass balances (m w.e.) due to the unknown density composition of snow, firn and ice.

September. All surface mass balances are provided in water equivalent rates (m w.e. y−1) averaged over the entire glacier basins. More information about the mass balance records of Kongsvegen, Hansbreen and Austfonna can be found in Hagen et al. (1999), Jania & Hagen (1996) and Moholdt et al. (2010), respectively. 4. Methods We tested three methods (Fig. 4) for estimating elevation changes from ICESat data: (a) crossover points (Section 4.1), (b) DEM-projected repeat-tracks (Section 4.2), and (c) planes fitted to repeat-tracks (Section 4.3). Section 4.4 describes how point data of elevation change were extrapolated to entire glacier regions in order to estimate regional volume changes and area-averaged elevation changes. Error analysis and validation of the methods follow in Section 5, while the glaciological results are presented and discussed in Section 6. For each of the three methods, we calculated average 2003–2008 elevation change rates (dh/dt), as well as seasonal winter and summer elevation changes (dhw and dhs) between the observation campaigns in fall (Oct./Nov.) and winter (Feb./Mar.). Note that the elevation change “winter season” (Oct./Nov.–Feb./Mar.) and “summer season” (Feb./ Mar.–Oct./Nov.) differ from the meteorological mass balance seasons which have a longer winter accumulation period and a shorter summer

A crossover point is the intersection between an ascending track and a descending track (Fig. 4a). Elevations at crossover points were linearly interpolated from the two closest footprints within 200 m in each track (e.g. Brenner et al., 2007). We used elevation differences at crossover points in three ways: (1) to estimate the RMS precision (σcross) of ICESat crossover data by only comparing crossovers within the same observation campaign (dt b 35 days) where only small elevation changes are expected, (2) to estimate average 2003–2008 elevation change rates (dh/dt) by comparing crossovers from similar seasons with a temporal separation of 3 or 4 years (dt≥ 3 y), and (3) to estimate seasonal elevation changes between fall (Oct./Nov.) and winter (Feb./Mar.) observation campaigns. Seasonal elevation changes (dhw and dhs) were estimated for each winter and summer season within the 2003–2008 ICESat data set. The crossover-point results were mainly used to validate the two repeat-track methods which provide a denser sample of elevation change points at the cost of a lower accuracy. 4.2. Elevation changes along DEM-projected repeat-tracks The unmeasured topography between near repeat-tracks needs to be considered when comparing elevations from different tracks. Slobbe et al. (2008) used a DEM to correct for the surface slope between the center points of overlapping footprints on the Greenland ice sheet. Using only overlapping footprints limits the slope-induced error, but it also limits the amount of data available for comparison. We applied a method which uses along-track interpolation to restrict the DEM slopecorrection to the cross-track distance between two repeat-tracks (Moholdt et al., 2010). For pairs of repeat-tracks, one profile is projected onto the other profile using the corresponding cross-track elevation differences from an independent DEM (Fig. 4b). Elevations are then compared at each DEM-projected point by linear interpolation between the two closest footprints in the other profile. The average cross-track separation between pairs of repeat-tracks at Svalbard was 73 m after removing repeat-track pairs separated by more than 200 m. Average annual elevation change rates (dh/dt) were calculated from all profile pairs obtained in similar seasons (i.e. winter–winter, summer–summer or fall–fall) with a temporal separation of 2, 3 or 4 years (dt≥ 2 y). Moholdt et al. (2010) used a minimum time span of 3 years to minimize short-term meteorological variations in the dh/dt estimates. In this study, we also included 2-year time spans to expand the spatial coverage in cloudy regions like Southern Spitsbergen where the amount of ICESat data is limited. Hence, most reference tracks include dh/dt points covering several different time spans. All these

Fig. 4. Three methods used to calculate elevation changes from ICESat data: (a) linear interpolation of neighbour footprints to crossover points between ascending and descending tracks (dh = HA–HB), (b) cross-track DEM projection (HDREF = HD2 + dHDEM) and linear interpolation to compare two repeat-tracks (dh = HDREF–HCREF), and (c) fitting least-squares regression planes to repeat-track observations to estimate slopes and average dh/dt.

2760

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

points were averaged within clusters every 350 m along-track to obtain mean 2003–2008 dh/dt values at a homogeneous spatial resolution. Winter and summer elevation changes (dhw and dhs) were calculated between the Oct./Nov and Feb./Mar. observation campaigns for the five ICESat winters from 2003/2004 to 2007/2008 and the four ICESat summers from 2004 to 2008. Along-track clustering was applied to each season.

plane residual and time for the earlier season (e.g. October 2006). The first term of Eq. (2) denotes the predicted elevation change under the assumption of a constant elevation change rate (dh/dt), while the second term adds the residual elevation change due to seasonal fluctuations in dynamics and accumulation/ablation.

4.3. Elevation changes at planes fitted to repeat-tracks

All elevation change points are aligned along a limited number of ICESat reference tracks (Fig. 1). In order to estimate ice volume changes, we need to extrapolate the observations to the remaining glacier areas. The data coverage of ICESat alone is typically too sparse over glaciers in mountainous terrain to allow local spatial interpolation. Therefore, we used a hypsometric approach to extrapolate over regions large enough to ensure that the distribution of ICESat tracks relative to individual glaciers is random, reducing the risk of systematic errors resulting from the measurement locations (Nuth et al., 2010). The relationship between elevation and dh/dt was parameterized for each region by fitting polynomial functions to the data (e.g. Kääb, 2008). The r2 coefficient of determination and the RMS error of the polynomial fits were typically stabilizing after adding a third order coefficient. Thus, we used third order polynomial fits p(h) to the elevation changes in each region:

Ideally, a DEM or other external data should not be required to compare near repeat-track elevations. A set of repeat-tracks contain a mixed elevation signal from local topography and temporal elevation changes between the observations. Several methods have been proposed to separate elevation changes from topographic variations using ICESat data only. Fricker & Padman (2006) computed elevation anomalies relative to averaged reference profiles, while Pritchard et al. (2009) compared elevation points to triangles spanned by three observations obtained within 2 years. We used a least-squares regression technique that fits rectangular planes to segments of repeat-track ICESat data (Howat et al., 2008). Along each reference track, multi-temporal ICESat points were assigned to 700 m long planes (Fig. 4c) overlapping by 350 m (i.e. most ICESat points belong to two planes). The width of the planes depends on the maximum cross-track separation distance between the repeated profiles, typically a few hundred meters. For each plane we estimated east and north slopes (αE, αN) and a constant elevation change rate (dh/dt) from the least-squares solution of the equation: 2

3 2 dE1 dH1 4 ⋮ 5=4 ⋮ dHn dEn

dN1 ⋮ dNn

3 32 2 3 r1 αE dt1 ⋮ 5⋅4 αN 5 + 4 ⋮ 5 rn dh=dt dtn

ð1Þ

where dE, dN, dH and dt are the differences in position and time (in decimal years) between each point and the average of all points in the plane. The residuals (r) of the plane regression contain remaining elevation variations which can not be ascribed to the assumption of planar slopes and an invariable elevation change rate. To avoid gross errors in dh/dt due to cloud-affected signals or small-scale topography, we removed potential outlier points where r N 5 m and recomputed the regression iteratively until all residuals were below this threshold. Finally, we removed all planes that consisted of less than 4 repeattracks or less than 10 points, as well as planes with a shorter observational time span than 2 years. The average number of tracks and points per plane after the filtering was 6.3 and 22, respectively. The time span of each plane varies greatly due to the scattered spatial coverage of each observation campaign. We anticipated that dh/dt estimates from planes with a different start and end season would be biased. For example, a plane spanning the whole repeat-track period from October 2003 to March 2008 would be slightly biased towards a more positive elevation change rate since the time span contains one more winter season than summer season. Therefore, the ICESat points in each plane were first filtered such that each plane start and end with the same season, i.e. winter to winter, summer to summer, or fall to fall. The filter was designed to obtain the longest possible time span with the highest number of points in the plane. We also investigated winter and summer elevation changes (dhw and dhs) from the plane regression results. For each of the five winter seasons (Oct./Nov.–Feb./Mar.) and four summer seasons (Feb./Mar.– Oct./Nov.) between fall 2003 and winter 2008, the seasonal elevation change (dhseas) of a plane was estimated from: P

P

P

P

dhseas = ð t seas1 − t seas0 Þ⋅dh=dt + ð r seas1 − r seas0 Þ

ð2Þ

where r̅seas1 and t̅seas1 are the average plane residual and time for the later season (e.g. March 2007), and r̅seas0 and t̅seas0 are the average

4.4. Volume changes and mass balance

3

2

pðhÞ = a1 ⋅h + a2 ⋅h + a3 ⋅h + a4

ð3Þ

where the a parameters are the least-squares solution to the function. Polynomial curves were fitted regionally to the winter elevation changes (dhw), the summer elevation changes (dhs) and the overall elevation change rates (dh/dt), with separate fits for the two repeattrack methods (Fig. 5). We also calculated polynomial fits for the 9 individual 2003–2008 winter and summer seasons, though it was necessary to combine all Svalbard data to obtain a sufficient sample. Volume changes (dV) were calculated from the equation: z

dV = ∑ ðpðhz Þ⋅Az Þ

ð4Þ

1

where p(h) is the third order polynomial function (Eq. (3)) fitted to the elevation changes, hZ is the middle elevation of 50 m elevation bins (e.g. 75 m for the 50–100 m elevation bin), and AZ is the glacier area for each of the Z elevation bins (Fig. 5). The 50 m glacier hypsometries were extracted from the glacier DEMs originating from between 1966 and 2008 depending on region. Similar volume changes were obtained if the mean (Arendt et al., 2002) or the median (Abdalati et al., 2004) of each bin was used instead of the polynomial fits. Area-averaged elevation changes were simply estimated by dividing the volume changes (dV) by the corresponding glacier areas (A). We calculated area-averaged winter and summer elevation Z Z changes (dhw and dhs) as well as area-averaged annual elevation change rates (dh=dt) (Table 1). We did not consider the influence of glacier area changes on the area-averaged elevation changes (e.g. Arendt et al., 2002) since we lack information about area changes within the 2003–2008 ICESat period. The regions in Fig. 1 do not include the glaciers of the Nordenskiöld peninsula (NP) in central Spitsbergen (720 km2) and a few smaller ice caps in the surroundings of Vestfonna and Austfonna (600 km2). In order to incorporate these sparsely sampled glaciers and ice caps in the overall Svalbard change estimates, we estimated the volume Z Z changes of the Nordenskiöld glaciers using the mean dhw , dhs and dh =dt of Northwestern and Southern Spitsbergen (Nuth et al., 2007), and the volume change of the smaller ice caps in the northeast from Z Z the mean dhw , dhs and dh =dt of Vestfonna and Austfonna. These extrapolated volume changes were then added to the total volume P

P

P

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

2761

Fig. 5. Third order polynomial fits to the elevation change rates (dh/dt) from the DEM method (red line) and the plane method (blue line) for the 7 glacier regions (Fig. 1) and for the entire Svalbard. The Svalbard subplot also includes a polynomial curve for the crossover dh/dt points (green line and dots). Dashed lines indicate the dh/dt variation within 50 m elevation bins represented as one standard deviation lines from the mean dh/dt of bins with at least 15 observations. Grey bars show the glacier hypsometries as area per 50 m elevation bin in the glacier DEMs. The lowermost lines represent the number of dh/dt observations per elevation bin for the DEM method (red) and the plane method (blue).

change of the 7 major regions to obtain the total Svalbard volume change, and the overall area-averaged elevation change (Table 1). Geodetic mass balances are typically estimated by multiplying dV/ dt with the ~0.9 density ratio between ice and water under the assumption of a constant firn pack. This simplification has usually little impact on decadal mass balance estimates (e.g. Nuth et al., 2010), but over short time spans, the effect of firn pack changes can be significant (e.g. Moholdt et al., 2010). In order to provide an estimate and uncertainty range for the overall 2003–2008 geodetic mass balance of Svalbard, we applied three simple density conversion schemes to the regional elevation change curves (Fig. 5): (1) assuming Sorge's Law (Bader, 1954) of a constant firn pack through time such that all changes are multiplied by the density of ice (ρice = 900 kg m−3), (2) assuming that all thinning consists of ice (ρice) while all thickening consists of firn (ρfirn ~ 500 kg m−3), and (3) assuming the density of ice (ρice) for changes in the lowermost 1/3 of the elevation bins (mainly ablation area) and the density of firn (ρfirn) for the uppermost 1/3 of the elevation bins (mainly accumulation area), with a linearly decreasing density from ρice to ρfirn for the middle 1/3 of the elevation bins. 5. Error analysis and validation of methods There is little external 2003–2008 elevation data to compare with, so the main validation was done by comparing the three elevation change methods with each other. The accuracy and precision of ICESat

elevation data have been thoroughly documented in other studies (e.g. Shuman et al., 2006), and crossover-point analysis has proved to be a very accurate way to compare ICESat elevations (e.g. Brenner et al., 2007). We used crossover points to estimate the error of individual repeat-track elevation change estimates, and then we used the RMS error of the polynomial fits to asses the error budget of the regional area-averaged elevation changes. The elevation change rates (dh/dt) from the three methods span a range of different time windows of 2–4 integer years (±35 days). They are thus influenced by mass balance variations during the 2003– 2008 period as well as short term snow variability within the ~35 days observation campaigns. The temporal data distribution (Fig. 2) does not point towards any particular observation campaigns that are heavily under- or over-represented. All in all, the many random time spans help to smooth out anomalous variations within the survey period. 5.1. Errors in the crossover points The RMS error of crossover-point comparisons (σcross) was estimated to be 0.66 m by comparing 329 glacier crossover points within individual ICESat observation campaigns (dt b 35 days) where only small elevation changes are expected. Outliers were removed through a 3σ filter which was run iteratively until the improvement of σcross was less than 5%. The error of crossover-point elevation change rates (σdh/dt) averaged over 3- or 4-year time spans is thus 0.22 or

2762

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

Table 1 Z Svalbard regions and their associated glacier surface area and volume change rate dV/dt. Area-averaged elevation changes are given for the 2004–2008 winter seasons dhw (Oct./ Z Nov.–Feb./Mar.), the 2003–2007 summer seasons dhs (Feb./Mar.–Oct./Nov.), and the 2003–2008 average annual dh =dt . All results are obtained using the repeat-track plane method. Glacier region Northwestern Spitsbergen (NW) Northeastern Spitsbergen (NE) Southern Spitsbergen (SS) Barentsøya and Edgeøya (BE) Vestfonna ice cap (VF) Austfonna ice cap (AF) Kvitøyjøkulen ice cap (KV) Regions total (REG) Svalbard total (SVAL)

Z

Z

Area (km2)

dhw (m)

dhs (m)

6300 8630 4760 2680 2410 7800 700 33,280 34,560

0.62 ± 0.15 0.61 ± 0.10 0.82 ± 0.19 0.87 ± 0.21 0.30 ± 0.13 0.46 ± 0.07 0.42 ± 0.25 0.60 ± 0.05 0.60 ± 0.05

−1.09 ± 0.16 −0.50 ± 0.10 −0.87 ± 0.23 −1.02 ± 0.29 −0.55 ± 0.17 −0.40 ± 0.08 −1.01 ± 0.30 −0.70 ± 0.06 −0.70 ± 0.06

dh =dt (m y−1)

dV/dt (km3 y−1)

−0.54 ± 0.10 0.06 ± 0.06 −0.15 ± 0.16 −0.17 ± 0.11 −0.16 ± 0.08 0.11 ± 0.04 −0.46 ± 0.11 −0.12 ± 0.04 −0.12 ± 0.04

−3.40 ± 0.72 0.52 ± 0.52 −0.71 ± 0.76 −0.46 ± 0.30 −0.39 ± 0.20 0.86 ± 0.32 −0.32 ± 0.08 −3.90 ± 1.27 −4.17 ± 1.28

The associated geodetic mass balance (excluding calving front fluctuations) of Svalbard is estimated to be −4.3 ± 1.4 Gt y−1, corresponding to an area-averaged balance of −0.12 ± 0.4 m w.e. y−1.

0.17 m y−1, with an average σdh/dt of 0.20 m y−1 (Table 2). The crossover error is generally lower in gentle terrain and higher in steeper or rougher terrain, with σcross progressively increasing from 0.34 m at 0–1˚ slopes (91 crossovers) to 0.86 m at 3–5˚ slopes (47 crossovers). 5.2. Errors in the DEM projections When a DEM is used to correct repeat-track ICESat data for the cross-track slope, it is not the absolute accuracy of the DEM that is important, but rather the reproduction of the relative local topography that is used to correct for the cross-track slope. We estimated the relative height error of each glacier DEM by calculating how well a DEM can predict the along-track elevation difference (i.e. the local slope) between neighbouring ICESat points separated by ~ 170 m. Since each ICESat elevation is obtained from a ~ 70 m diameter footprint, we anticipated that DEM smoothing would improve the correspondence with alongtrack ICESat slopes. For each DEM, we applied an iterative low pass mean filter of increasing pixel size (i.e. 3 × 3, 5 × 5, 7 × 7 etc.) until the improvement of the RMS of the ICESat–DEM point-pair differences was less than 5%. The optimal averaging window sizes for the Spitsbergen SPOT DEMs were 7 × 7 (280 × 280 m) in NW and SS, and 9 × 9 (360 × 360 m) in Northeastern Spitsbergen. The DEMs interpolated from contour maps and the InSAR/ICESat DEM at Austfonna were already sufficiently smooth with no significant improvement from additional averaging. The final RMS values ranged from 1.0 m at Austfonna to 2.5–3.5 m in the semi-alpine Spitsbergen regions after applying an iterative 3σ filter with a convergence threshold of 5%. These values are upper estimates of the DEM-projection error (σDEM) since the cross-track separations are typically much less than 170 m. In addition to σDEM comes the along-track interpolation error which should be less than the crossover-point error (σcross). The combined dh/dt error varies greatly in space depending on the repeattrack separation distance (0–200 m), the quality of the DEM, the length of the time span, as well as the surface slope and roughness. The error is reduced through the averaging of elevation change points within 350 m clusters. We estimated the error of DEM-projected elevation changes rates (σdh/dt) by comparing the crossover dh/dt points (dt ≥ 3 y) with the closest of these clusters within a 500 m radius. This resulted in 307 comparable dh/dt points with an RMS error of 0.48 m y−1 after applying the iterative 3σ filter (Fig. 6). The errors of the seasonal elevation change points from the winter (σwinter) and summer (σsummer) periods were similarly estimated to be 1.09 m and 1.21 m from samples of 193 and 130 points, respectively (Table 2).

5.3. Errors in the plane fitting The application of the plane method assumes that the regression scheme (Eq. (1)) is able to separate between the slopes of a plane (αE and αN) and the average elevation change rate (dh/dt). The along-track slope component (α||) is typically well resolved by each repeat-track, while the cross-track slope component (αX) of a plane is dependent on a number of noncoincident repeat-tracks which are influenced by dh/dt. Comparisons between crossover points and their closest plane within a 500 m radius show that the two data sets mostly yield consistent estimates of dh/dt (Fig. 6) and αX (Fig. 7) with no signs of systematic errors. The elevation change error (σdh/dt) of the plane method was estimated to be 0.34 m y−1 from the RMS of the dh/dt differences at 294 crossovers which remained after the iterative 3σ filter (Fig. 6). Similarly, the errors of the winter (σwinter) and summer (σsummer) elevation changes were estimated to be 0.78 and 0.93 m from samples of 194 and 130 crossover points, respectively (Table 2). We also compared the plane cross-track slopes (αX) with the corresponding slopes extracted from the surrounding 3 × 3 pixels in the smoothed version of the DEMs. This allowed us to validate

Fig. 6. Validation of repeat-track elevation change estimates (dh/dt) close to crossoverpoint locations. About 300 crossover dh/dt points (dt ≥ 3 y) are compared to the closest repeat-track dh/dt point within 500 m distance for the DEM method and the plane method. The RMS errors yield the estimated dh/dt accuracies (σdh/dt) for the two repeat-track methods (Table 2).

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767 Table 2 Estimated RMS errors (σ) for individual estimates of elevation change. All uncertainties are based on data comparisons at a few hundred crossover-point locations at Svalbard glaciers. Method

σwinter

σsummer

σdh/dt

Crossovers DEM method Plane method

0.66 m 1.09 m 0.78 m

0.66 m 1.21 m 0.93 m

0.20 m y−1 0.48 m y−1 0.354 m y−1

al errors for the average 2003–2008 elevation change rates (ε̄dh/dt) and for the average seasonal elevation changes during winter (ε̄w) and summer (ε̄s) are shown in Table 1. Regional volumetric errors (E) were obtained from the root-sumsquares (RSS) of the specific error (ε̄) multiplied with the regional glacier area (A) and the volume change (dV) multiplied with a tentative glacier area uncertainty of ±10% (Berthier et al., 2010): E=

all ~ 9000 planes against an independent data set. The RMS differences of αX were 1.24˚ for the DEM comparison and 0.58˚ for the crossover-point comparison (Fig. 7). Although the noise level of the DEMs is high, the results confirm that the plane method is able to resolve the local topography in a reasonable way. 5.4. Errors in the area-averaged elevation changes Area-averaged elevation change errors (ε ) were estimated for each region and method by applying the standard error equation on each polynomial fit:  ε

σfit = pffiffiffiffiffiffiffiffiffiffiffi N−4

ð5Þ

where σfit is the RMS error of the polynomial fit, N-4 is the degrees of freedom in the third order polynomial function with 4 unknown parameters (Eq. (3)), and N is the number of uncorrelated elevation change observations in the region. It is not straightforward to quantify correlation distances and magnitudes for elevation change measurements (e.g. Rolstad et al., 2009). Most error assessments are based on simplified assumptions about the spatial autocorrelation. Nuth et al. (2010) pointed out that individual ICESat profiles are correlated, but restricted the correlation distance to within 50 m elevation bins. Moholdt et al. (2010) assumed that all elevation change observations on Austfonna were fully correlated within 2 km clusters, while the clusters themselves were uncorrelated. This study has fewer altimetric data sources and a larger potential for measurement correlation. Therefore, we chose a conservative correlation distance of 5 km (~ 14 planes or clusters) for the dh/dt observations. Thus, N in Eq. 5 refers to the number of 5 km along-track segments containing elevation change data. The estimated region-

2763

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð ε ⋅AÞ2 + ðdV⋅0:1Þ2

ð6Þ

The RSS of the regional volumetric errors (E) form the overall Svalbard volumetric error, assuming no data correlation between the regions. The area-averaged elevation change error for the entire Svalbard is then found by dividing the overall volume change error by the total glacier area. In this way, we estimated the uncertainty of the area-averaged 2003–2008 elevation change rate (ε̄dh/dt) at Svalbard to ±0.04 m y−1 for both repeat-track methods (Table 1). When a limited sample of elevation change points is used to estimate regional glacier changes, we can divide the error budget into an observation error (ε̄OBS) and a spatial extrapolation error (ε̄EXT) (e.g. Arendt et al., 2002; Nuth et al., 2010). We modified Eq. (5) to estimate the area-averaged observation error (ε̄OBS): σ  ffiffiffiffiffiffiffiffiffiffiffi ε OBS = pmethod N−4

ð7Þ

where σmethod represents the methodological elevation change errors in Table 2, and N-4 is the degrees of freedom for N uncorrelated 5 km segments of elevation change data. The resulting ε̄OBS estimates for the Z plane method at Svalbard are ±0.01 m y−1 for dh=dt, ±0.04 m for dhw Z and ±0.03 m for dhw , with slightly higher uncertainties for the DEM method. The magnitude of ε̄OBS will vary from region to region depending on the length of ICESat profiles and the along-track topography, but there are too few crossover points to estimate one σmethod for each region. However, the relative height errors of the DEMs (Section 5.2) indicate that ε̄OBS is about twice as high for the semi-alpine Spitsbergen regions as for the more gentle ice caps on the other islands. We assume that the ICESat profiles are randomly distributed in space such that the spatial elevation change variations are captured by the observations and thus included in the overall area-averaged errors. Knowing the overall error (ε̄) and the observation error (ε̄OBS), we can then estimate the area-averaged extrapolation error (ε̄EXT) from the RSS relation between them:  ε EXT =

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ε 2 − ε 2OBS

ð8Þ

The extrapolation errors (ε ̅EXT) Zat Svalbard are ±0.03 m y−1 Z (dh=dt ) and ±0.04 m (dhw and dhs ) for the plane method, and slightly lower for the DEM method which have more dh/dt data than the plane method (Fig. 5). In order to check if the errors are realistic, we randomly assigned the plane tracks into two independent data sets and calculated separate change rates and errors. The two sets of results were within the error bounds of each other for all regions, and the Svalbard dh= dt values were within ±0.01 m of the original estimate using all data. A case study at Edgeøya also confirms that a limited number of ICESat profiles can yield a good estimate of the overall glacier change as compared to multi-temporal DEMs (Kääb, 2008).

P

5.5. Comparison of the three elevation change methods

Fig. 7. Validation of cross-track slope estimates (αx) of planes. The plane αx values are compared to corresponding slopes of the closest crossover point within 500 m distance and to slopes calculated from the smoothed glacier DEMs at the plane locations.

The crossover-point method (Fig. 4a) provides the most accurate elevation change points (Table 2). The overall polynomial elevation change curve from crossover points is similar to those of the two repeattrack methods (Fig. 5), but the spatial coverage is too sparse for regional

2764

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767 P

volume change calculations. We found that if the Svalbard dh=dt was calculated from one of the three overall Svalbard elevation change curves, then it would be ~0.05 m y−1 higher than the sum of the regional changes in Table 1. This bias is due to a spatial under-sampling of the thinning regions in the west and south where cloud cover and rugged topography sometimes hinder the elevation change calculations. At the regional scale, we assume that the spatial sampling of dh/dt is random and thus has no impact on the final results. The comparison of the dh/dt rates of the two repeat-track methods with neighbouring crossover points, shows good agreement despite a high noise level (Fig. 6). The error estimates of the plane method are lower than for the DEM method (Table 2), implying that the topographic signal within ICESat repeat-tracks is actually more precise than the existing DEMs at Svalbard. The polynomial dh/dt curves from the two methods agree well in all regions, with no signs of systematic differences (Fig. 5). The one-standard-deviation curves (dashed lines) are generally closer to the dh/dt curves for the plane method than for the DEM method. The larger per-point errors in the dh/dt estimates of the DEM method are compensated by a better spatial coverage, yielding a similar area-averaged error ( ε̅dh/dt) as the plane method. A potential problem with the plane method for dh/dt calculation is the uneven temporal data sampling (Fig. 2). There is typically more data from the winter campaigns, and less data from the summer campaigns. The risk of a seasonal bias in dh/dt is especially high for planes where the earliest and latest ICESat observations stem from different seasons. If all available ICESat data were to be included in the plane calculations, the regional polynomial curves would be shifted upwards by 0.05– 0.10 m y−1 as compared to the curves of the DEM method, which should not be biased since it only compares data between similar seasons. This bias towards more positive dh/dt results is probably due to the fact that the overall ICESat epoch from October 2003 to March 2008 spans one more winter season than summer season. After applying the seasonal plane filter (Section 4.3), the overall Svalbard dh=dt fell from −0.06 to −0.12 m y−1 which is similar to the DEM method. All the final results in Table 1 are based on the plane method which provides elevation change data purely from ICESat. 6. Results and discussion 6.1. Multi-year elevation changes Fig. 1 shows the average annual 2003–2008 elevation change rates (dh/dt) at plane and cluster locations. There are large local dh/dt variations which we attribute to differences in glacier dynamics, wind drift and measurement noise. Data in the western parts of Svalbard are generally noisier and more diverse than the eastern parts, making it more difficult to interpret elevation change trends in the west. All regions apart from Vestfonna are characterized by a low-elevation thinning of 1–3 m y−1 (Fig. 5). Elevation change rates at the higher elevations are generally more positive, though with variations from slight thinning in Northwestern Spitsbergen to pronounced thickening of up to 0.5 m y−1 in Northeastern Spitsbergen and at Austfonna. The general tendency of surface steepening implies that most Svalbard glaciers are not in dynamic equilibrium with the recent surface mass balance (Hagen et al., 2005). The strong frontal thinning has caused a significant glacier area shrinkage, a trend that has been observed over most of Svalbard (e.g. Nuth et al., 2007; Dowdeswell et al., 2008; Blaszczyk et al., 2009). Such geometric change patterns are typical for glaciers in the quiescent phase of their surge cycle (e.g. Melvold and Hagen, 1998). However, some of the surface steepening may be related to short term meteorological factors like firn area expansion and thickening which have probably occurred on Svalbard when the surface mass balance turned from very negative in 2003/ 2004 to more balanced between 2004 and 2007 (Fig. 3). The most striking feature in Fig. 1 is the extensive interior thickening in Northeastern Spitsbergen and at the Austfonna ice cap.

The recent thickening at Austfonna has been confirmed by repeated GNSS surface profiles (Moholdt et al., 2010), and the pattern is also consistent with 1996–2002 airborne laser profiles (Bamber et al., 2004). In Northeastern Spitsbergen, substantial thickening has previously only been observed in the southern part, more specifically in the accumulation area of the Negribreen basin which surged around 1935 (Hagen et al., 1993). The other high-elevation areas in the region have typically been in balance or been slightly thinning during 1965–2005 (Nuth et al., 2010) and 1996–2002 (Bamber et al., 2005). Both these studies concluded that the Northeastern Spitsbergen glaciers are generally less negative than the glaciers in Northwestern and Southern Spitsbergen. The glaciers in Northwestern Spitsbergen have recently been thinning more than in Southern Spitsbergen (Figs. 1 and 5), a trend which is also evident in the surface mass balance curves of Kongsvegen (NW) and Hansbreen (SS) (Fig. 3). This north–south tendency is opposite of what has been found from earlier geodetic data (Bamber et al., 2005; Nuth et al., 2007; Kohler et al., 2007; Nuth et al., 2010). The high elevations of the Southern Spitsbergen glaciers have recently been close to balance (Fig. 5), while thinning was dominating between 1990 and 2005 (Nuth et al., 2010). A similar shift towards less negative dh/dt seems to have occurred on the glaciers of Barentsøya and Edgeøya as compared to 1970–2002 elevation change curves from DEM differencing at the Kvalpyntfonna and Digerfonna ice caps (Kääb, 2008). The most homogenous elevation change pattern is seen on Austfonna, while the adjacent Vestfonna ice cap is the region with the most complex changes (Fig. 1). At Vestfonna, both thinning and thickening occur at all elevations depending on location. For example, the main summit has thinned by up to 0.5 m y−1, while the tidewater fronts of Franklinbreen to the northwest have thickened, and also advanced (Sneed, 2007). These change patterns have also been recognized between 1990 and 2005 (Nuth et al., 2010) and more recently from GNSS surface profiles and field observations (V. Pohjola, unpublished data). Most of the basin–scale complexity is likely caused by the glacier dynamics with surge-type and active surging outlet glaciers (Dowdeswell & Collin, 1990). Wind redistribution of snow might also have an impact (Beaudon & Moore, 2009). The local anomalous changes make the elevation change curve of Vestfonna to deviate from the other regions (Fig. 5). The Kvitøyjøkulen ice cap in the far east of Svalbard has generally thinned, especially at the lower elevations (Figs. 1 and 5). There are no published records of elevation changes at Kvitøyjøkulen, but crossover points between 1983 airborne radio echo-sounding (RES) profiles (Bamber & Dowdeswell, 1990) and ICESat profiles indicate an overall thinning over the last few decades. In addition, large volumes of ice have been lost along the ~100 km calving front which has retreated by an average of ~25 m y−1 between the 1983 survey and the recent ICESat measurements. 6.2. Seasonal elevation changes Seasonal elevation changes were calculated for the ICESat winter seasons (Oct./Nov.–Feb./Mar.) and summer seasons (Feb./Mar.–Oct./ Nov.) between October 2003 and March 2008. The errors of individual seasonal elevation change estimates are higher than for the multiyear elevation change rates due to the shorter time spans involved (Table 1). The overall area-averaged seasonal elevation changes for Z Svalbard are 0.60 ± 0.05 m for the average winter season ( dhw ) and Z −0.70 ± 0.06 m for the average summer season dhs ). The consistency of the elevation changes in Table 1 can be tested by comparing Z Z the sums of dhw and dhs with the corresponding annual elevation change rates (dh=dt). All differences are within the error bounds. The inter-regional differences in winter and summer elevation Z Z changes (dhw and dhs ) are often smaller than the associated uncertainties (Table 1), but a few spatial patterns can be recognized. The western and southern regions have a larger mass turnover than the P

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

northeastern regions, characterized by more thickening during winter (0.62–0.87 vs. 0.30–0.61 m) and more thinning during summer (0.87– 1.09 vs. 0.40–0.55 m). This pattern agrees well with long term regional mass balance gradients (Hagen et al., 2003b) as well as the more recent surface mass balance time series which are fluctuating more for Kongsvegen (NW) and Hansbreen (SS) in western Spitsbergen than for Etonbreen at Austfonna (Fig. 3). Note that the ICESat winter and summer seasons do not correspond to the meteorological accumulation and ablation seasons, implying that the real seasonal elevation changes should have larger magnitudes than those presented here. A cumulative time series of area-averaged Svalbard elevation changes for the 5 winter seasons and 4 summer seasons between 2003 and 2008 is shown in Fig. 3. The first year of record from October 2003 to October 2004 is responsible for an area-averaged thinning of roughly 0.6 m which is comparable to the total ice loss during 2003– 2007. Both the smallest winter thickening and the largest summer thinning occurred during the anomalous 2003/2004. The following 3 years from 2004 to 2007 are more balanced with comparable winter increases and summer decreases. Similar temporal trends can be recognized in the surface mass balance curves of Kongsvegen (NW), Hansbreen (SS) and Etonbreen (AF), although there are local variations (Fig. 3). The overall cumulative elevation change is slightly positive if all seasons are included and slightly negative if one of the extra winter seasons are kept out. The large seasonal and regional variations in elevation change and mass balance make it difficult to interpret the 2003–2008 changes in a longer time perspective. 6.3. Volume changes and mass balance The average annual volume change (dV/dt) of Svalbard glaciers during 2003–2008 is estimated to be −4.17 ± 1.28 km3 y−1 (excluding calving front fluctuations), corresponding to an area-averaged elevation change rate of −0.12 ± 0.04 m y−1 (Table 1). Most of the regions have experienced volume losses, with the largest change in Northwestern Spitsbergen at −3.40 ± 0.72 km3 y−1. Northeastern Spitsbergen and Austfonna are the only regions which have gained volume, totally 1.38 ± 0.61 km3 y−1. If we exclude Austfonna and Kvitøyjøkulen, the remaining Svalbard glaciers have thinned at an average rate of −0.14 ± 0.05 m y−1 which is much less than the −0.39 ± 0.02 m y−1 loss rate between 1965–1990 DEMs and 2003– 2008 ICESat (Nuth et al., 2010). Northeastern Spitsbergen has turned from an average thinning of −0.27 ± 0.02 m y−1 (1965–2005) to an average thickening of 0.06 ± 0.06 m y−1 (2003–2008), while the average thinning in Southern Spitsbergen and at Barents-/Edgeøya has slowed down from −0.55 ± 0.06 to −0.17 ± 0.11 m y−1. In contrast to this, Northwestern Spitsbergen and Vestfonna appear to be more negative recently, though the differences are barely significant due to the high uncertainties in these regions (Table 1). Note that the numbers above represent glacier volume changes and area-averaged elevation changes which can deviate by different amounts from the associated mass balances due to variations in the amount and density of snow and firn within the survey periods. The area-averaged 0.11 ± 0.04 m y−1 thickening at Austfonna is higher than what Moholdt et al. (2010) found by combining multitemporal 2002–2008 elevation data from GNSS surface profiles, airborne laser altimetry and ICESat. Most of the 0.06 m y−1 difference can be explained by the update of glacier outlines which removed 200 km2 of low-lying ice. The surface mass balance at Austfonna went from strongly negative in 2003/2004 to strongly positive in 2007/ 2008 (Moholdt et al., 2010) with the result of an expanding and thickening firn layer as observed by ground-penetrating radar (Dunse et al., 2009). Satellite ASTER imagery over Spitsbergen glaciers shows that the firn extent was much smaller in summer 2003 than in the summer 2007/2008 SPOT 5 scenes. This development has caused an increase in the amount of firn at Svalbard between 2003 and 2008, despite an overall glacier thinning. Therefore, regional mass

2765

balances are probably lower than they would be if the 0.9 ice/ water density ratio was used to convert the geodetic changes in Table 1 into mass balances. The overall geodetic mass balance of Svalbard was assessed by applying three simple density conversion schemes to the regional elevation change curves (Section 4.4). Conversion scheme (1) resulted in an overall mass balance of −3.7 Gt y−1, while scheme (2) and (3) yielded mass balances of −4.9 and −4.2 Gt y−1 respectively. Considering that the firn area has expanded between 2003 and 2008 and that firn layers densify with time (e.g. Reeh, 2008), we expect that scheme (1) and (2) are extreme cases of density conversion in either end. Hence, our estimated glacier volume change of −4.17 km3 y−1 (Table 1) will most likely correspond to a geodetic mass balance within the range [−3.7 – 4.9] Gt y−1. If we use the mean value of the range and combine the uncertainty with the dV/dt error as a RSS, it results in a 2003–2008 Svalbard geodetic mass balance of −4.3 ± 1.4 Gt y−1 or −0.12 ± 0.4 m w.e. y−1 when averaged over the total glacier area. Potential glacier area changes are not accounted for although glacier retreat has been observed in Southern Spitsbergen (Nuth et al., 2007), at Austfonna (Dowdeswell et al., 2008) and for most tidewater glaciers (Blaszczyk et al., 2009). As an effect of general area shrinkage, the volume change rates are probably slightly too negative for regions with glacier outlines prior to the ICESat epoch (BE, VF, and KV). The areaaveraged rates, on the other hand, are mostly influenced by glacier area changes within the survey period. This error is negligible for 2003–2008 if we assume an average annual area shrinkage of 0.3 % which was observed in Southern Spitsbergen between 1936 and 1990 (Nuth et al., 2007). However, glacier area changes at the terminus of an advancing or retreating tidewater glacier will also include ice volume changes below sea level, which cannot be measured by ICESat. Blaszczyk et al. (2009) used 2000–2006 ASTER imagery to estimate a total Svalbard calving front length of 860 ± 15 km and an average terminus retreat of 30 ± 10 m y−1 yielding a marine loss rate of −2.3± 0.8 Gt y−1 when assuming an average terminus ice thickness of 100 ± 10 m. Accounting for this, the total Svalbard mass balance during 2003–2008 becomes −6.6± 2.6 Gt y−1, corresponding to an area-averaged balance of −0.19 ± 0.08 m w.e. y−1. Previous Svalbard mass balance estimates from ~30 years of insitu measurements have been hampered by a lack of ablation measurements outside of western Spitsbergen. Hagen et al. (2003a; 2003b) found that the estimated surface mass balance could diverge between −10 and −0.5 Gt y−1 depending on whether western ablation data were used in all regions or whether regional mass balance curves were estimated from the available accumulation data, respectively. In addition comes the total ice loss from calving which has been estimated to be −4 ± 1 Gt y−1 (Hagen et al., 2003b) and more recently to −6.8 ± 1.7 Gt y−1 (Blaszczyk et al., 2009). Most combinations of these numbers yield a more negative mass balance than our recent estimate. The large regional variations detected by ICESat exemplify the difficulty of extrapolating surface mass balance data from a few locations to the entire Svalbard ice mass. The uncertainty of the ice flux calving rate is also lower in the geodetic results since glacier dynamics is an integrated part of the derived volume change. However, satellite altimeters are unable to measure density changes which can have a significant influence on the mass balance, especially over short time spans like in this study. The Svalbard contribution to sea-level change is not straight forward to quantify due to the uncertain distinction between marine ice losses above and below sea level. If we assume an average terminus ice thickness of 100 m (Hagen et al., 2003b) and an average ice cliff height of 25 m above sea level (derived from ICESat), then the contribution from terminus retreat or advance to sea-level change will be ~ 13 % of the water equivalent ice volume change at the terminus. Thus, Svalbard's −2.3 ± 0.8 Gt y−1 (Blaszczyk et al., 2009) terminus retreat loss will only displace 0.3 ± 0.1 Gt y−1 of sea-water,

2766

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767

yielding a 2003–2008 Svalbard contribution to sea-level change of 4.6 ± 1.4 Gt y−1, corresponding to 0.013 ± 0.004 mm sea-level equivalent. This is lower than a previous estimate of 8.8 ± 3 Gt y−1 from GRACE gravity measurements between February 2003 and January 2008 (Wouters et al., 2008). 7. Conclusions We have tested two ICESat repeat-track methods for calculating glacier elevation changes at the Svalbard archipelago. The derived elevation changes are sparse and not always reliable on a local scale, but hypsometric averaging within larger glacier regions provides robust estimates of volume changes and area-averaged elevation changes. The regional estimates benefit from various track configurations over many glaciers rather than center-line profiles over a few selected glaciers (e.g. Abdalati et al., 2004; Bamber et al., 2005) which can cause a bias due to spatial sampling issues (Berthier et al., 2010). The multi-temporal sampling of ICESat help to smooth out anomalous meteorological seasons which can have a large impact on the overall elevation change rates and mass balances over short time periods like 2003–2008. The DEM method and the plane method (Fig. 4) yield consistent results (Fig. 5) and agree well with more accurate elevation change calculations at crossover-point locations (Fig. 6). The good performance of the plane method implies that it can also be used in other Arctic regions of similar characteristics, even when glacier DEMs are not available. We recommend applying a seasonal data filter (Section 4.3) prior to the plane fitting if the goal is to determine average annual elevation change rates. This is to correct for the bias towards more positive elevation change rates due to the one additional winter season within the ICESat epoch. Most glacier regions in Svalbard have experienced low-elevation thinning and high-elevation balance or thickening (Fig. 5). The overall 2003–2008 geodetic mass balance (excluding calving front retreat or advance) is estimated to be −4.3 ± 1.4 Gt y−1, corresponding to an area-averaged thinning of −0.12 ± 0.4 m w.e. y−1. This is less negative than the previous few decades as estimated from comparing ICESat with 1965–1990 DEMs (Nuth et al., 2010). The largest areaaveraged thinning has happened in the west and south, while Northeastern Spitsbergen and the Austfonna ice cap have slightly thickened. The average summer thinning has been about twice as high in the west and south as in the northeast (Table 1), probably related to shorter ablation seasons in the northeast than in the west and south as observed by satellite scatterometer data (Rotschky et al., 2008). Time series of seasonal elevation changes can be obtained if the sparse seasonal ICESat observations are averaged over a sufficiently large glacier area. The cumulative seasonal elevation change curve for Svalbard shows that most ice losses occurred during the very negative 2003/2004 mass balance year (Fig. 3). The temporal trends derived from ICESat fit well with the surface mass balance curves of Kongsvegen (NW), Hansbreen (SS) and Etonbreen (AF). The large spatial and temporal variations in elevation change at Svalbard make it difficult to interpret the results from this study in a long term perspective. Future altimetry records from CryoSat-2 and ICESat-2 will be crucial for determining if the current trends are representative for the climatic development at Svalbard. Altimetric time series like these can also help to validate coincidental gravity measurements from the GRACE satellite which have been applied with success in the Gulf of Alaska glacier region (Luthcke et al., 2008; Arendt et al., 2008). Finally, there is a need to relate the 2003–2008 elevation changes with meteorological data as a tool for mass balance modelling both back in time and into the future. Acknowledgements The authors are very thankful to the numerous data contributors that made this study possible, namely the National Snow and Ice Data Center (ICESat data), the SPIRIT project (SPOT 5 stereoscopic survey of Polar Ice:

Reference Images and Topographies), the Norwegian Polar Institute (topographic maps and Kongsvegen mass balance), the Polish Academy of Sciences, P. Glowacki and D. Puczko (Hansbreen mass balance), the University of Oslo / NPI (Austfonna mass balance), and J. Dowdeswell and T. Benham (1983 RES data). Funding for this study was partly provided by the CryoSat calibration and validation experiment (CryoVEX) coordinated by the European Space Agency, and the International Polar Year project GLACIODYN; the dynamic response of Arctic glaciers to global warming. The final stage was supported by funding to the ice2sea project from the European Union 7th Framework Programme, grant number 226375, ice2sea contribution number 009. G. Moholdt was also supported through the Arktisstipend grant from the Svalbard Science Forum (SSF). Furthermore, we acknowledge A. Kääb, B. E. Smith and two anonymous reviewers for useful comments and suggestions which helped to improve the manuscript.

References Abdalati, W., Krabill, W., Frederick, E., Manizade, S., Martin, C., Sonntag, J., Swift, R., Thomas, R., Yungel, J., & Koerner, R. (2004). Elevation changes of ice caps in the Canadian Arctic Archipelago. Journal of Geophysical Research-Earth Surface, 109. Abshire, J. B., Sun, X. L., Riris, H., Sirota, J. M., McGarry, J. F., Palm, S., Yi, D. H., & Liiva, P. (2005). Geoscience Laser Altimeter System (GLAS) on the ICESat mission: On-orbit measurement performance. Geophysical Research Letters, 32. Arendt, A., Luthcke, S. B., Larsen, C. F., Abdalati, W., Krabill, W., & Beedle, M. J. (2008). Validation of high-resolution GRACE mascon estimates of glacier mass changes in the St Elias Mountains, Alaska, USA, using aircraft laser altimetry. Journal of Glaciology, 54, 778−787. Arendt, A. A., Echelmeyer, K. A., Harrison, W. D., Lingle, C. S., & Valentine, V. B. (2002). Rapid wastage of Alaska glaciers and their contribution to rising sea level. Science, 297, 382−386. Bader, H. (1954). Sorge's Law of densification of snow on high polar glaciers. Journal of Glaciology, 2, 319−323. Bamber, J., Krabill, W., Raper, V., & Dowdeswell, J. (2004). Anomalous recent growth of part of a large Arctic ice cap: Austfonna, Svalbard. Geophysical Research Letters, 31. Bamber, J. L., & Dowdeswell, J. A. (1990). Remote-sensing studies of Kvitøyjøkulen, an ice cap on Kvitøya, North-east Svalbard. Journal of Glaciology, 36, 75−81. Bamber, J. L., Krabill, W., Raper, V., Dowdeswell, J. A., & Oerlemans, J. (2005). Elevation changes measured on Svalbard glaciers and ice caps from airborne laser data. Annals of Glaciology, 42, 202−208. Beaudon, E., & Moore, J. (2009). Frost flower chemical signature in winter snow on Vestfonna ice cap, Nordaustlandet, Svalbard. The Cryosphere, 3, 147−154. Berthier, E., & Toutin, T. (2008). SPOT5-HRS digital elevation models and the monitoring of glacier elevation changes in North-West Canada and South-East Alaska. Remote Sensing of Environment, 112, 2443−2454. Berthier, E., Schiefer, E., Clarke, G. K. C., Menounos, B., & Remy, F. (2010). Contribution of Alaskan glaciers to sea-level rise derived from satellite imagery. Nature Geoscience, 3, 92−95. Bjornsson, H., Gjessing, Y., Hamran, S. E., Hagen, J. O., Liestol, O., Palsson, F., & Erlingsson, B. (1996). The thermal regime of sub-polar glaciers mapped by multi-frequency radio-echo sounding. Journal of Glaciology, 42, 23−32. Blaszczyk, M., Jania, J. A., & Hagen, J. O. (2009). Tidewater glaciers of Svalbard: recent changes and estimates of calving fluxes. Polish Polar Research, 30, 85−142. Brenner, A. C., DiMarzio, J. R., & Zwally, H. J. (2007). Precision and accuracy of satellite radar and laser altimeter data over the continental ice sheets. Ieee Transactions on Geoscience and Remote Sensing, 45, 321−331. Dowdeswell, J. A., & Collin, R. L. (1990). Fast-flowing outlet glaciers on Svalbard ice caps. Geology, 18, 778−781. Dowdeswell, J. A., Benham, T. J., Strozzi, T., & Hagen, J. O. (2008). Iceberg calving flux and mass balance of the Austfonna ice cap on Nordaustlandet, Svalbard. Journal of Geophysical Research-Earth Surface, 113. Dunse, T., Schuler, T. V., Hagen, J. O., Eiken, T., Brandt, O., & Høgda, K. A. (2009). Recent fluctuations in the extent of the firn area of Austfonna, Svalbard, inferred from GPR. Annals of Glaciology, 50, 155−162. Fricker, H. A., Borsa, A., Minster, B., Carabajal, C., Quinn, K., & Bills, B. (2005). Assessment of ICESat performance at the Salar de Uyuni, Bolivia. Geophysical Research Letters, 32. Fricker, H. A., & Padman, L. (2006). Ice shelf grounding zone structure from ICESat laser altimetry. Geophysical Research Letters, 33. Fricker, H. A., Scambos, T., Bindschadler, R., & Padman, L. (2007). An active subglacial water system in West Antarctica mapped from space. Science, 315, 1544−1548. Førland, J., & Hanssen-Bauer, I. (2003). Past and future climate variations in the Norwegian Arctic: overview and novel analyses. Polar Research, 22, 113−124. Hagen, J. O., Liestøl, O., Roland, E., & Jørgensen, T. (1993). Glacier atlas of Svalbard and Jan Mayen. Meddelser nr. 129. Oslo: Norwegian Polar Institute 141 pp. Hagen, J. O., Melvold, K., Eiken, T., Isaksson, E., & Lefauconnier, B. (1999). Mass balance methods on Kongsvegen, Svalbard. Geografiska Annaler Series a-Physical Geography, 81A, 593−601. Hagen, J. O., Kohler, J., Melvold, K., & Winther, J. G. (2003). Glaciers in Svalbard: mass balance, runoff and freshwater flux. Polar Research, 22, 145−159.

G. Moholdt et al. / Remote Sensing of Environment 114 (2010) 2756–2767 Hagen, J. O., Melvold, K., Pinglot, F., & Dowdeswell, J. A. (2003). On the net mass balance of the glaciers and ice caps in Svalbard, Norwegian Arctic. Arctic, Antarctic, and Alpine Research, 35, 264−270. Hagen, J. O., Eiken, T., Kohler, J., & Melvold, K. (2005). Geometry changes on Svalbard glaciers: mass-balance or dynamic response? Annals of Glaciology, 42, 255−261. Hamilton, G. S., & Dowdeswell, J. A. (1996). Controls on glacier surging in Svalbard. Journal of Glaciology, 42, 157−168. Howat, I. M., Smith, B. E., Joughin, I., & Scambos, T. A. (2008). Rates of southeast Greenland ice volume loss from combined ICESat and ASTER observations. Geophysical Research Letters, 35. Hutchinson, M. F. (1989). A new procedure for gridding elevation and stream line data with automatic removal of spurious pits. Journal of Hydrology, 106, 211−232. Jania, J., & Hagen, J. O. (1996). Mass Balance of Arctic Glaciers. Report No. 5. SosnowiecOslo: International Arctic Science Committee 62 pp. Johannessen, O. M., Khvorostovsky, K., Miles, M. W., & Bobylev, L. P. (2005). Recent icesheet growth in the interior of Greenland. Science, 310, 1013−1016. Kaser, G., Cogley, J. G., Dyurgerov, M. B., Meier, M. F., & Ohmura, A. (2006). Mass balance of glaciers and ice caps: consensus estimates for 1961–2004. Geophysical Research Letters, 33. Kohler, J., James, T. D., Murray, T., Nuth, C., Brandt, O., Barrand, N. E., Aas, H. F., & Luckman, A. (2007). Acceleration in thinning rate on western Svalbard glaciers. Geophysical Research Letters, 34. Korona, J., Berthier, E., Bernard, M., Remy, F., & Thouvenot, E. (2009). SPIRIT. SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies during the fourth International Polar Year (2007–2009). Isprs Journal of Photogrammetry and Remote Sensing, 64, 204−212. Kääb, A. (2005). Remote Sensing of Mountain Glaciers and Permafrost Creep. Schriftenreihe Physische Geographie (266). : University of Zürich. Kääb, A. (2008). Glacier Volume Changes Using ASTER Satellite Stereo and ICESat GLAS Laser Altimetry. A Test Study on Edgeoya, Eastern Svalbard. Ieee Transactions on Geoscience and Remote Sensing, 46, 2823−2830. Lefauconnier, B., & Hagen, J. O. (1991). Surging and calving glaciers in Eastern Svalbard. Meddelelser. Oslo: Norwegian Polar Institute 130 pp. Luthcke, S. B., Arendt, A., Rowlands, D. D., McCarthy, J. J., & Larsen, C. F. (2008). Recent glacier mass changes in the Gulf of Alaska region from GRACE mascon solutions. Journal of Glaciology, 54, 767−777. Melvold, K., & Hagen, J. O. (1998). Evolution of a surge-type glacier in its quiescent phase: Kongsvegen, Spitsbergen, 1964–95. Journal of Glaciology, 44, 394−404. Moholdt, G., Hagen, J. O., Eiken, T., & Schuler, T. V. (2010). Geometric changes and mass balance of the Austfonna ice cap, Svalbard. The Cryosphere, 4, 21−34. Muskett, R. R., Lingle, C. S., Sauber, J. A., Rabus, B. T., & Tangborn, W. V. (2008). Acceleration of surface lowering on the tidewater glaciers of Icy Bay, Alaska, USA from InSAR DEMs and ICESat altimetry. Earth and Planetary Science Letters, 265, 345−359. Nuth, C., Kohler, J., Aas, H. F., Brandt, O., & Hagen, J. O. (2007). Glacier geometry and elevation changes on Svalbard (1936–90): a baseline dataset. Annals of Glaciology, 46, 106−116. Nuth, C., Moholdt, G., Kohler, J., Hagen, J. O., & Kääb, A. (2010). Svalbard glacier elevation changes and contribution to sea level rise. Journal of Geophysical Research, 115. Pinglot, J. F., Hagen, J. O., Melvold, K., Eiken, T., & Vincent, C. (2001). A mean net accumulation pattern derived from radioactive layers and radar soundings on Austfonna, Nordaustlandet, Svalbard. Journal of Glaciology, 47, 555−566. Pritchard, H. D., Arthern, R. J., Vaughan, D. G., & Edwards, L. A. (2009). Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets. Nature, 461, 971−975.

2767

Reeh, N. (2008). A nonsteady-state firn-densification model for the percolation zone of a glacier. Journal of Geophysical Research-Earth Surface, 113. Rolstad, C., Haug, T., & Denby, B. (2009). Spatially-integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: application to the Western Svartisen ice cap, Norway. Journal of Glaciology, 55, 666−680. Rotschky, G., Kohler, J., & Isaksson, E. (2008). Spatio-temporal variability of snow melt on Svalbard derived from spaceborne scatterometer data (Qscat). Workshop on The Dynamics and Mass Budget of Arctic Glaciers, 29–31 January 2008, Obergurgl, Austria. Extended abstracts. IMAU Utrecht: IASC Working Group on Arctic Glaciology. Sauber, J., Molnia, B., Carabajal, C., Luthcke, S., & Muskett, R. (2005). Ice elevations and surface change on the Malaspina Glacier, Alaska. Geophysical Research Letters, 32. Schutz, B. E., Zwally, H. J., Shuman, C. A., Hancock, D., & DiMarzio, J. P. (2005). Overview of the ICESat Mission. Geophysical Research Letters, 32, 4. Shuman, C. A., Zwally, H. J., Schutz, B. E., Brenner, A. C., DiMarzio, J. P., Suchdeo, V. P., & Fricker, H. A. (2006). Shuman, C.A., Zwally, H.J., Schutz, B.E., Brenner, A.C., DiMarzio, J.P., Suchdeo, V.P., & Fricker, H.A. (2006). ICESat Antarctic elevation data: Preliminary precision and accuracy assessment. Geophysical Research Letters, 33. Geophysical Research Letters, 33. Slobbe, D. C., Lindenbergh, R. C., & Ditmar, P. (2008). Estimation of volume change rates of Greenland's ice sheet from ICESat data using overlapping footprints. Remote Sensing of Environment, 112, 4204−4213. Smith, B. E., Bentley, C. R., & Raymond, C. F. (2005). Recent elevation changes on the ice streams and ridges of the Ross Embayment from ICESat crossovers. Geophysical Research Letters, 32. Smith, B. E., Fricker, H. A., Joughin, I. R., & Tulaczyk, S. (2009). An inventory of active subglacial lakes in Antarctica detected by ICESat (2003–2008). Journal of Glaciology, 55, 573−595. Sneed, W. (2007). Satellite Remote Sensing of Arctic Glacier - Climate Interations. M. Sc. thesis, 83 pp. University of Maine. Sund, M., Eiken, T., Hagen, J. O., & Kääb, A. (2009). Svalbard surge dynamics derived from geometric changes. Annals of Glaciology, 50, 50−60. Thomas, R., Frederick, E., Krabill, W., Manizade, S., Martin, C., & Mason, A. (2005). Elevation changes on the Greenland ice sheet from comparison of aircraft and ICESat laser-altimeter data. Annals of Glaciology, 42, 77−82. Wingham, D. J., Ridout, A. J., Scharroo, R., Arthern, R. J., & Shum, C. K. (1998). Antarctic elevation change from 1992 to 1996. Science, 282, 456−458. Wingham, D. J., Francis, C. R., Baker, S., Bouzinac, C., Brockley, D., Cullen, R., de ChateauThierry, P., Laxon, S. W., Mallow, U., Mavrocordatos, C., Phalippou, L., Ratier, G., Rey, L., Rostan, F., Viau, P., & Wallis, D. W. (2006). CryoSat: a mission to determine the fluctuations in Earth's land and marine ice fields. Natural Hazards and Oceanographic Processes from Satellite Data, 37, 841−871. Wise, S. (2000). Assessing the quality for hydrological applications of digital elevation models derived from contours. Hydrological Processes, 14, 1909−1929. Wouters, B., Chambers, D., & Schrama, E. J. O. (2008). GRACE observes small-scale mass loss in Greenland. Geophysical Research Letters, 35. Zwally, H. J., Brenner, A. C., Major, J. A., Bindschadler, R. A., & Marsh, J. G. (1989). Growth of the Greenland ice sheet - measurement. Science, 246, 1587−1589. Zwally, H. J., Schutz, B., Abdalati, W., Abshire, J., Bentley, C., Brenner, A., Bufton, J., Dezio, J., Hancock, D., Harding, D., Herring, T., Minster, B., Quinn, K., Palm, S., Spinhirne, J., & Thomas, R. (2002). ICESat's laser measurements of polar ice, atmosphere, ocean, and land. Journal of Geodynamics, 34, 405−445. Zwally, H. J., Schutz, R., Bentley, C., Bufton, J., Herring, T., Minster, B., Spinhirne, J., & Thomas, R. (2008). GLAS/ICESat L1B Global Elevation Data V028, 20 February 2003 to 21 March 2008. Boulder. CO: National Snow and Ice Data Center Digital media.