Decadal changes from a multi-temporal glacier inventory of Svalbard

Introduction Conclusions References


Introduction
Glacier inventories are important for studying the global frozen freshwater resource, and provide a basic dataset for further glaciological, remote sensing and modeling studies.The World Glacier Inventory (WGI), the first global glacier catalogue, was compiled with classification schemes based on hydrological drainage basins (Müller et al., 1977).WGI contains auxiliary information such as topographic parameters, length and volume estimates, and glacier characterization codes, but it does not include the digitized coordinates of the glacier outlines (WGMS, 1989).Recently, the Global Land Ice Measurements from Space (GLIMS) initiative has provided a scheme for generating a global glacier inventory that retains the raw glacier outline information (Raup et al., 2007).There are some inherent differences between GLIMS and WGI in their structure, application and information they provide (Cogley, 2009), but regional glacier inventories can be relatively easily submitted and linked to both datasets (e.g.Paul et al., 2009;Bolch et al., 2010;Svoboda and Paul, 2010).Moreover, the increasing ease of generating glacier inventories semi-automatically from satellite imagery (Paul et al., 2013) combined with the more frequent coverage of planned satellite missions (e.g.Sentinel-2) facilitate such investigations of multi-temporal glacier inventory changes.
The Arctic archipelago of Svalbard in the North Atlantic is ∼ 57 % glacierized and contains a mix of cirque and valley glaciers, ice fields and ice caps.There are both landterminating and tidewater glaciers; most of them are polythermal and many exhibit surge-type behavior.The surface mass balance has generally been negative since the termination of the Little Ice Age, which ended around Svalbard in the beginning of the 1920s; by this time most glaciers had reached their maximum Neocene extent (Hagen et al., 1993(Hagen et al., , 2003)).Summer temperatures increased dramatically during the 1920s and 1930s (Nordli and Kohler, 2004) in this part

C. Nuth et al.: A multi-temporal glacier inventory of Svalbard
of the Arctic, a period sometimes referred to as the Early 20th Century Warming (Wood and Overland, 2010).Following a colder period from the 1940s to the 1960s, summer temperatures on Svalbard have been increasing.For the period 1989-2011, summer temperatures have been increasing at rates of more than 0.5 • C decade −1 at the two long-term meteorological stations (Førland et al., 2011).This has led to increased glacier volume loss, particularly in western Svalbard (e.g.Kohler et al., 2007;James et al., 2012).
In this study, we present glacier extent snapshots and change rates from the 1930s to 2010 based on 3 digital glacier inventories: GI old , GI 90 and GI 00s .The inventories are a key component of a new digital glacier database (König et al., 2013) that is structured after the Glacier Atlas of Svalbard and Jan Mayen (Hagen et al., 1993), the first complete glacier inventory of the archipelago (referred to as H93 in the rest of this text).GI old and GI 90 are derived from two Norwegian Polar Institute map products; the first is a mixed product of the years 1936, 1960, 1961, 1966, 1969, 1970, and 1971, and the second is from 1990.GI 00s updates the previous inventories using satellite imagery from the years 2000-2010.We describe the present glaciation of the archipelago through topographic and glaciologic inventory parameters.Furthermore, we discuss the generation and applicability of three glacier inventory change parameters as derived for two epochs: (1) area changes, (2) length changes as estimated by two methods, and (3) glacier tongue width change.

Historic data
Accurate topographic mapping of Svalbard began in the 1950s, when the Norwegian Polar Institute (NPI) undertook to construct maps using analogue photogrammetry of oblique aerial photographs taken in 1936 and 1938.These early maps covered southern and western Spitsbergen, about 22 % of the archipelago.From the late 1950s to the early 1970s, a number of aerial campaigns acquired vertical photographs covering ∼ 50 % of the archipelago, with major campaigns in 1966 for northern Spitsbergen and in 1971 for Edgeøya/Barentsøya.Taken together, these maps formed the original S100 (1 : 100 000) topographic map series published and distributed by NPI as paper maps.The original S100 map series was digitized by NPI in the 1990s, and forms the basis for GI old .
H93, the original glacier inventory of Svalbard (Hagen et al., 1993), followed the identification and parameter definitions outlined by WGMS (1989).It was based upon the S100 paper maps (before digital transformation) but with the oldest data (1936 and 1938) updated using pre-1980 aerial and Landsat imagery.Front positions and areas of H93 were measured from these updated paper maps using a planimeter.The inventory consists of basic data such as glacier name, area, length, and photo year in table format, but the raw outline locations are not preserved and thus not available digitally in a GIS (geographic information system).For consistency, we preserve the same structure for our new digital glacier inventory.
A second major mapping campaign was conducted by NPI in 1990 to acquire vertical aerial photographs over nearly the entire archipelago.In the 1990s and 2000s, NPI created topographic and thematic maps for about 60-70 % of the archipelago using digital photogrammetric techniques and manual feature delineation, which forms the basis for GI 90 .There are two exceptions in this updated dataset: the south coast of Austfonna front position was mapped by helicopter with GPS in 1992 (T. Eiken, personal communication, 2013), and a small inland strip within southern Spitsbergen was covered with 1961and 1970images (H. Faste Aas, personal communication, 2013).The spatial and temporal coverage of all glacier inventories can be seen in Fig. 1.

Satellite imagery
Summer satellite imagery spanning the period 2000-2010 is used as the basis for our third glacier dataset, GI 00s (Table 1).We prioritize data from sensors that obtain stereo optical imagery for creation of orthophotos that are temporally and spatially consistent with the digital elevation models (DEMs) used to generate them.Accordingly, the main data are the DEMs and orthophotos from the SPOT5-HRS (highresolution sensor) satellite, generated by the IPY project SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies) (Korona et al., 2009).The SPOT5-HRS collects 5 m panchromatic stereo images that are stereoscopically processed into 40 m DEMs, which are then used for generating orthophotos of the original images.Five SPIRIT acquisitions from 2007 to 2008 cover 71 % of the glacierized area of Svalbard.
The second main satellite data source is the ASTER L1B product, in the form of automatically generated DEMs and orthophotos (AST14DMO, 2010).These have a smaller swath width (60 km) than SPOT5, such that 23 scenes are needed to cover ∼ 16 % of the glacier area.Cloud-free scenes were not available for 2007-2008, and therefore data from as early as 2000 were required to complete coverage of the archipelago.For the remaining 14 % of the glacier area, no suitable SPOT5-HRS or ASTER scenes were available.For these glaciers, 11 orthorectified Landsat scenes from 2001 to 2007 are used.An additional 17 Landsat and 13 ASTER scenes are used to aid decisions related to perennial snow patches and to delineate glacier margins in areas with seasonal snow cover or shade.

Satellite DEMs
We use the ASTER GDEM (v2) for separating glaciers into individual hydrological units and for prescribing topographic attributes to the glaciers (Frey and Paul, 2012).The GDEM is a global compilation of stacked and filtered ASTER DEMs (Fujisada et al., 2012) that temporally overlaps GI 00s .GDEM is chosen over other DEMs (SPOT/ASTER) for its consistency as a single product covering the entire archipelago.Glacier surfaces in the GDEM have a bumpy texture, a result of the merging of temporally different DEMs of varying quality, especially on the low visual contrast upper glacier areas.Therefore, a low-pass Fourier filter is applied over glacier surfaces to remove the high-frequency noise and minimize the size of the blunders that occur at the highest elevations (see Appendix A).The filtering reduced standard deviations of elevation differences with respect to ICESat and SPOT-SPIRIT DEMs and improved visual appearance of the GDEM without changing the overall structure of the surface.Moreover, visual comparisons between the GDEM-derived hydrological basins and those derived from the NPI topographic maps/DEMs, the SPOT-SPIRIT DEMs and individual AST14DMO (2010) DEMs reveal small variations which verify the use of the GDEM for this purpose and infer that rough DEM quality does not have a large impact on drainage basin generation.The largest discrepancies (blunders) occur on the flattest upper regions of ice cap and ice fields where small elevation inconsistencies can sometimes lead to large differences in the determination of a hydrological divide.These blunders are manually adjusted to the drainage basins derived from the other DEMs.For the Austfonna ice cap, where the GDEM contains several holes, we use an independently created DEM (Moholdt and Kääb, 2012), as well as velocity fields derived from SAR interferograms (Strozzi et al., 2008;Dowdeswell et al., 2008) to delineate ice divides and glacier basins and to generate topographic parameters.

Georeferencing
The various DEMs and orthoimages must be correctly georeferenced for merging into a common dataset.While the SPOT5 and ASTER orthoimages are internally consistent with the associated DEMs, the geolocation accuracy is dependent upon the accuracy of the satellite position determination (orbital parameters) and instrument pointing (auxillery attitude information), and thus the relative DEM/orthophoto may not necessarily be located precisely on the ground.We co-register the DEMs to ICESat laser altimetry data (Zwally et al., 2012) over non-glacier topography (Nuth and Kääb, 2011) and apply the horizontal component of the correction vector to the orthoimages.We reference ICESat rather than the NPI S100 map series because ICESat data are acquired in a consistent way over the entire archipelago, while S100 is a merged product originating from various sources (analogue and digital photogrammetry) and dates (1960s-1990s).
The horizontal consistency of ICESat to the 1990 DEM is ≈ ±3 m (Nuth and Kääb, 2011).For the LANDSAT scenes, manual co-registration is performed by applying a linear translation based on tie points from NPI coastline vector data and the available co-registered satellite imagery.

Glacier delineation and identification
The raw outlines that form the historic inventories of GI old and GI 90 were generated by cartographers who visually interpreted and digitized the glacier borders from the original images using analogue photogrammetric workstations for the oldest dataset (GI old ) and digital orthoimages for GI 90 .These outlines were highly detailed (high resolution) with accurate glacier front positions; however, many of the raw outlines contained seasonal snow-covered valley walls and gullies higher up on the glaciers from misinterpretation by the cartographers.These were trimmed by visual analysis of the multi-temporal satellite archives where it was obvious that the cartographers digitized snow-filled gullies, which are not considered a part of the glacier and were removed by best judgment rather than using a minimum size criteria.These historic glacier outlines did not distinguish between individual glaciers, but rather were complex polygons.We combine the raw digitized glacier outlines from S100 maps  with the analogue H93 inventory to create GI old .Glacier basins are delineated based on the H93 local identification codes (WGI IDs) and glacier names, and with the visual help of all available automatically generated hydrological basins, topographic contours and optical satellite imagery (Sect.2).GI old is then used as a reference to separate glacier basins (Racoviteanu et al., 2009) in the 1990 outlines, forming the GI 90 inventory.Finally, the GI 00s dataset is created from the most recent of either GI 90 or GI old by manually trimming or reshaping the front position and the lateral edges of the glacier tongue to the more recent satellite orthophotos.The basis of the outline reshaping is visual interpretation of the satellite optical images using contrast stretching to aid the delineation of the glacier margin.Outlines were also updated in the upper regions of the glaciers when nunataks appeared or large changes were present due to upper glacier downwasting, for example from a surge.The local identification system for individual glaciers is defined hydrologically by the WGI IDs of H93 comprising 5 digits: -1st digit represents the division of the archipelago into 5 major regions: (1) Spitsbergen, (2) Nordaustlandet, (3) Barentsøya, (4) Egdeøya, and (5) Kvitøya.
-2nd digit is the division of each region into primary drainage basins.
-3rd digit is the division into secondary drainage basins.
-4th and 5th digits are the number for each individual glacier.
For example, if a glacier is denoted by 14 204, then the glacier lies in region 1 (Spitsbergen), in major drainage basin 4 (Isfjorden), and in secondary drainage basin 2 (Adventdalen), and its glacier number is 04 (Longyearbreen).An overview map of the regions and drainage basins is shown in Fig. 2. The original H93 glacier identification system required adaptation since a number of individual WGI glacier units in H93 comprised single tongues fed by multiple tributary ice streams that can now be divided into new discrete flow units.Either the old glacier front has retreated and naturally separated into separate tongues or the tributary glaciers have significant medial moraines suggesting a natural division of the glacier system.These new independent glacier units retain the original H93 ID, but with additional decimals to identify the individual glacier entities (see Fig. 3).In addition, a few of the flow divides on the 2 larger ice caps have been adjusted based upon updated information (Sect.2.3) that was not available in the creation of H93.
Lastly, an inventory requires also definitions of the smallest snow patches and glacierets (Cogley et al., 2011).GI 00s defines glacierets and snow patches as those visually and perennially present in the temporal series of Landsat and ASTER images, the smallest of which is 0.05 km 2 .H93 defines glacierets and snowpatches as perennial snow/ice areas less than 1 km 2 , and these units were not given specific IDs.Therefore, these units are labeled in GI 00s based on their secondary drainage basin, using 99 as the glacier number (4th and 5th digits) with increasing decimals for each individual unit (14 299.01 in the above example).13116   1 3 1 1 7 . 1   13118  13201   13202.1  1 3 2 0 2 . 2  1 3 2 0 2 . 3   13 20 5

Glaciological and topographic attributes
For the GI 00s inventory, a number of descriptive, glacier and topographic attributes are extracted for each individual glacier entity, as suggested by Paul et al. (2009).The standard geometric and topographic parameters include polygon area, polygon perimeter, elevation minimum, maximum, median and standard deviation, mean slope and mean aspect.Glacier hypsometries are extracted for 50 m elevation bins from the ASTER GDEM; these are shown in Fig. 2 for the primary drainage basins.Two additional glaciological parameters are generated for the three inventories: glacier length and the average width of the glacier tongue.At least one centerline is manually digitized for each glacier area polygon, from the glacier tongue to the head of the accumulation area.If a single glacier entity contains more than one centerline, the maximum length is provided.The centerlines are then used to generate glacier tongue width.Lines perpendicular to the centerline are intersected with the glacier outlines for each GI.The glacier tongue or terminus width is then estimated as the average measured width of the lowermost 10 % of the centerline for GI 00s .The threshold is chosen visually to best represent the varying tongue shapes of both small and large glaciers collectively.Varying the threshold by 5 % has little effect except for the smallest glaciers with pointy glacier tongue shapes.For GI 90 and GI old , if the centerline length change is greater than 10 % of the earlier centerline, the average width along the area of change is used to ensure estimates are representative for the area of change within an epoch.For glaciers containing multiple centerlines, the average front width of all centerlines are provided.For glaciers that have 2 separate tongues corresponding to individual centerlines, the sum of front widths is given.Glacier tongue widths serve two purposes in the inventory.The first is to estimate a calving front width and the second is for change analysis.Examples of the basic geometric structures of the inventories are shown in Fig. 3 and a list of attributes is given in Table 3.

Glacier change parameters
A common parameter for comparing multi-temporal inventories is area change (e.g.Le Bris et al., 2011;Davies and Glasser, 2012), expressed either as mean annual area change or relative change rates as a percentage.Analysis of raw and relative area changes alone is complicated by its dependence on glacier area, tongue width and other geometrical parameters.Therefore, we also derive length changes by two methods.First, centerline length change rates are calculated as the difference in centerline length between the inventories, divided by the time interval.Second, we use an adaptation of the "box method" employed in Greenland (Moon and Joughin, 2008;Box and Decker, 2011;Howat and Eddy, 2011), which provides an average change across the glacier tongue rather than a single estimate dependent upon the location of the centerline.In our approach, an average length change rate is defined as the area change below the GI 00s median elevation divided by the oldest glacier tongue width (as described above) and subsequently by the time interval.We refer to this length change estimate as the "area/width" length change.
Each GI contains glacier outlines from multiple times (Fig. 1).This complicates the interpretation of changes between the entire inventories.However, each individual glacier outline contains a single date associated with an image.Therefore, all change rates provided are specifically calculated for each glacier's individual time separation between outlines.Moreover, entire inventory comparisons are made by first calculating the change rates for individual outlines and then summing to drainage basin or region.

Inventory characteristics
The newest inventory, GI 00s , contains 1668 individually labeled glacier units (including snowpatch and glacier IDs with a decimal) totaling 33 775 km 2 , or about 57 % of the total land area of the archipelago.The distributions of glacier lengths and sizes are approximately log-normal; there are about 350 glaciers (22 % of the inventory) larger than 10 km 2 that make up 93 % of the glacier area (Fig. 4a).Similarly, there are 630 glaciers, glacierets and snowpatches smaller than 1 km 2 that represent 38 % of the inventory and 1 % of the glacier area (Table 2).Glacier centerline lengths range from 200 m to 60 km, with an average of 4.5 km.A significant loglinear (power-law) relationship exists between glacier area and length, as predicted by theory (Bahr, 1997) and shown with global inventory data (Bahr and Radić, 2012).
About 68 % of the glacierized area drains through tidewater calving fronts; their spatial distribution is shown in Fig. 6.The exact number of tidewater glaciers depends on how a glacier is defined.H93 labeled 152 glaciers with a calving front.In GI 00s , 197 glaciers with unique identification codes (including decimals) are characterized as tidewater terminating.Twenty-nine of these individually labeled glacier units share the calving glacier tongue with other glaciers but are clearly divided by medial moraines; this leaves 168 calving fronts.Blaszczyk et al. ( 2009) classified 163 tidewater glaciers, with the difference of 5 glaciers comprising 11 glaciers not classified as tidewater in GI 00s and 6 glaciers not classified as tidewater in Blaszczyk et al. (2009).
Summing the estimated front widths for tidewater glaciers results in a calving front length of 740 km, about half (386 km) of which are tidewater fronts in Spitsbergen.Our glacier tongue widths represent flux gates rather than the precise calving front length.Since lateral edges of many tidewater glaciers in Spitsbergen are grounded on land, our width estimates may often be larger than the active dynamical flux gate of the glacier.The difference with another calving front length estimate of 860 km (Blaszczyk et al., 2009) is due to our smaller front widths on the lobate tongues of ice cap outlet glaciers in Nordaustlandet (215 km), Edgeøya (23 km) and Kvitøya (113 km).
The islands to the east and northeast of Spitsbergen (Edgeøya, Barentsøya and Nordaustlandet) contain the flattest topographies, and thus glaciers there have lower slopes (Fig. 5b) and are mostly characterized by ice cap geometries.Consequently, glacier hypsometries typically feature an abrupt truncation at the highest elevations (Fig. 2).These islands contain about 40 % of the glacierized area of Svalbard, with heavy glaciation on Nordaustlandet (60-90 %) and to a lesser extent on Barentsøya/Edgeøya (Fig. 4b).Maximum and median elevations are lower for the ice caps of Barentsøya and Edgeøya (Figs. 2 and 5a).About 80 % of the Nordaustlandet glacier area drains through tidewater calving fronts, while only 47 % of the Barentsøya/Edgeøya glacier area is tidewater (Fig. 4d).
Spitsbergen contains more alpine topography than the islands to the east.Ahlmann et al. (1933) described the "Spitsbergen style" glacier as "continuous ice masses divided into individual ice streams by mountain ridges and nunatak areas".Spitsbergen glacier hypsometries are more normally distributed than those of the ice caps to the northeast (Fig. 2).The area-elevation distributions are positively skewed, indicating a greater hypsometric weight towards lower elevations.Three clusters of interconnected ice fields exist in northwest, northeast and south Spitsbergen.These ice field clusters are divided by a less glacierized central and northcentral interior (Fig. 4b).This central region contains the largest numbers of glaciers (Fig. 4c) with the highest median elevations and steeper average slopes (Fig. 5a, b), reflecting the more cirque-style glaciation in these alpine areas.The three ice field clusters contain all the tidewater glaciers of  Table 2. Glacier statistics for the major drainage basins of Svalbard for H93 and GI 00s (in bold).Total area is the size of the drainage basin, both glacier and land."Glacierized area" is the total glacierized area in each atlas.Differences in this parameter include glacier change and ommission/commission errors between the datasets."Comparable glacier area" is the area corresponding to similar IDs in both inventories excluding snow patches and glacierets (see Sect. 3.2).Differences in this parameter include glacier changes and an error associated with delineating glaciers."Percent area change" is derived from "omparable glacier area".The number of glacier units is provided only for GI 00s as the number of unique 5-digit IDs (no decimals) provides the total number for H93.For GI 00s , this is the number of merged integer IDs.Also shown is the number of individual snow patches (GI 00s IDs = XXX99.XX) and glaciers less than 1 km (H93) along with the area sums.All area estimates have unit km 2 .Spitsbergen (Fig. 6), which together drain about 62 % of the Spitsbergen glacierized area.The main tidewater drainage occurs off the eastern and western coastline of Spitsbergen and in Hornsund, Van Keulenfjord, Kongsfjord and Krossfjord.

Comparison of complete glacier inventories (H93 vs. GI 00s )
This section describes the differences between the only two fully complete glacier inventories of Svalbard, namely H93 and GI 00s .The comparison is complicated by both the varying time spans from which they were derived (Fig. 1) and the lack of raw H93 outlines to control that the upper glacier boundaries are consistent between the inventories.H93 contains 1022 individually labeled glaciers and 1164 non-labeled snow/ice masses less than 1 km 2 totaling 36 633 km 2 , or 60 % of the archipelago's land area (Table 2), while GI 00s contains 1668 individual glacier area polygons, or ∼ 57 % of the archipelago land area.GI 00s contain more individual units (polygons) due both to glacier retreat and separation, as well as our identification of distinct glacier flow units in glaciers previously classified by a single ID.Thus, GI 00s glaciers are combined to their parent single 5digit integer ID for comparisons.GI 00s contain an additional 52 smaller glaciers not formally identified in H93; these we have defined with 5-digit integer IDs that continue from the highest glacier number in each secondary drainage basin.Note that they do not follow the standard counter-clockwise identification sequence of H93.
In terms of snow patches, GI 00s contains roughly a quarter of the number and area of snow patches as H93.There existed a large number of thin snow polygons in the NPI historic maps.These are due to local topographic depressions (i.e.gullies, trenches, chutes) that often remain snow-filled in mid-to late summer, and thus were liable to be identified by cartographers as "glacier" in aerial and satellite images.We do not consider these polygons as glaciers and they are not maintained in our GI inventories; however, they are most likely included in the H93 count of glaciers/snow patches < 1 km 2 and probably are the reason behind the differences in snowpatches between the inventories.The difference between comparable areas of H93 and GI 00s reveals a glacierized area loss of 7 % (Table 2), or about 80 km 2 a −1 for the average 32 yr time span between the inventories.Area changes are computed at the primary and secondary drainage basin scale (Table 2 and Fig. 6b).This reduces random uncertainties from individual glacier divides and possibly misclassifications (see Sect. 4.4).The smallest relative changes (−2 to −5 %) have occurred in Nordaustlandet and the largest (−13 to −17 %) in central Spitsbergen, a region dominated by small glaciers, and Barentsøya/Edgeøya (Table 2).These patterns naturally reflect glacier area itself (Figs.4b and 6b) as described in many other glacier inventory studies (Kääb et al., 2002;Andreassen et al., 2008;Racoviteanu et al., 2008;Bolch et al., 2010;Le Bris et al., 2011) and further complicate spatial and temporal change analysis of inventory data.

Comparison of multi-temporal glacier inventories
(GI old vs. GI 90 vs. GI 00s ) The digital glacier inventories are available for three different times (Fig. 1), with 946 consistent glaciers (∼ 30 % of the total Svalbard glacierized area) located in southern and western Spitsbergen.This permits analysis of two time periods, from GI old to GI 90 (Epoch 1, which is ∼ 50 yr on average) and from GI 90 to GI 00s (Epoch 2, which is ∼ 17 yr on average).In sum, these glaciers lost ∼ 31 km 2 a −1 (0.26 % a −1 ) during Epoch 1 and ∼ 24 km 2 a −1 (0.23 % a −1 ) during Epoch 2.
In the following analysis, the sample population is limited to glaciers larger than 2 km 2 (406 glaciers) since smaller glaciers are more prone to interpretation errors related to seasonal snow.
Figure 7a shows centerline length changes and area/width length changes (calculated according to Sect.3.4).On average, area/width retreat rates are ∼ 10 m a −1 larger than centerline retreat rates.This is expected since the area/width length changes are an integrated change that also include lateral losses, while the centerline length changes are dependent upon one measurement taken along the centerline at the glacier front.
The distribution of length change rates (as estimated using the area/width approach) in both epochs is shown in  −150 m a −1 , with an average of ∼ −40 m a −1 and median of −30 m a −1 for both epochs.Larger extreme retreat rates exist (±150-350 m a −1 ), in most cases related to surge behavior.It seems the retreat of glaciers has increased, as is apparent in the shift in the distribution of retreat rates towards more negative values in Epoch 2. This is mostly apparent in central and southern Spitsbergen, as compared to northwest Spitsbergen, where average length change rates have remained similar (or even less negative) between the epochs (Fig. 8).
Comparisons between the relative changes of Epoch 1 and Epoch 2 for the area and length changes display varying patterns (Fig. 7b).Relative length changes are mainly larger for Epoch 2 than Epoch 1. Similarly to the absolute differences described above (Fig. 7a), relative centerline changes are smaller and vary more than relative area/width length changes.Relative area changes show greater scatter, with many glaciers experiencing larger relative changes during Epoch 1 than Epoch 2. This pattern is opposite that of relative length changes (both centerline and area/width length changes) and is at least partially a result of lateral glacier wastage, which was larger in Epoch 1 than Epoch 2. Nevertheless, all relative changes are dependent upon the original size of the parameter (length or area), and thus spatial and temporal comparisons are hampered by this dependence, which results in heteroscedastic distributions with glacier size (see e.g.Kääb et al., 2002;Bolch et al., 2010).

Accuracy
Errors in the glacier outlines depend on the images used to delineate glaciers (i.e.their resolution and quality) and sky and ground conditions, and the analyst's ability to digitize and interpret the imagery.Errors of the latter kind arise both from the manual interpretation of glacier-land boundaries and from the uncertainty of locating hydrological divides of interconnected ice fields (i.e. based upon surface topography).Errors in ice field divides are related to the accuracy of the DEM and to the hydrological flow directions derived from it when using automated hydrological GIS algorithms.Interpretation uncertainty may arise, for example, in cases where debris or lateral moraines obscure the glacier outlines, or where seasonal snow in the imagery covers the glacier edge.A manual digitization experiment (Paul et al., 2013) with 20 participants on 24 glaciers resulted in area uncertainties (expressed as a relative difference) ranging between 2 and 30 %; the largest errors came from sections of glaciers with heavy debris cover.Manual digitization error was found to be on the order of 1-3 pixels at any vertex; relative errors were typically better than 5 %, varying with glacier size and conditions (i.e.debris cover) (Paul et al., 2013).For our digital datasets, we expect errors of this magnitude but also some degree of spatial variability in the uncertainty since, for example, central and north-central areas are less glacierized (i.e. less than 40 % in Fig. 4b) and have larger amounts of debris cover and/or ice-cored moraines.
Glacier outlines for H93 are not digitally available, but they are based on many of the same topographic maps as GI old from which we have derived glacier divides independently using historic and recent DEMs (Sect.3).Glacier areas based on data from the same year can be compared to estimate an uncertainty related to glacier division and manual delineation.There are 170 common glacier units in H93 and GI old ; their relative differences approximate a Student's t distribution (i.e.heavier tails) with a standard deviation of about 8 %, while the standard deviation of the Gaussian distribution fit is 20 % (Fig. 9).The heavier tails of this relative error distribution result from gross differences in determining www.the-cryosphere.net/7/1603/2013/The Cryosphere, 7, 1603-1621, 2013 drainage divides or from the inclusion or exclusion of lateral moraines, which impacts the relative error more heavily than the smaller random errors introduced in digitization (as described in Paul et al., 2013).We define the individual glacier area error as the 95 % confidence interval of the Student t distribution, about 16 %, but note that the relative error is dependent upon glacier size as well (Kääb et al., 2002;Bolch et al., 2010;Paul et al., 2013), with smaller glaciers having larger relative errors.The bulk of the glaciers have errors less than 5 %, similar to other inventories (Paul et al., 2002;Bolch et al., 2010;Gjermundsen et al., 2011;Rastner et al., 2012).The error may be largely systematic at the individual glacier scale but is random at the regional or inventory-wide scale; i.e. the uncertainty of the drainage divides is cancelled.A rough conservative estimate for the error of the entire glacierized area of Svalbard is 1-2 % (∼ 500 km 2 ).
Finally, we simulate a 16 % error on the area changes and on the glacier tongue widths to estimate a sensitivity to and the precision of our area/width length changes.For the entire population of changes from Epoch 2, the residuals between the original length changes and those calculated with 16 % differences in the area changes and widths separately result in a error distribution approximated by a Student t with 95 % of the residuals contained within ±10 m a −1 .Combining both uncertainties from area changes and glacier tongue width es- Fig. 9. Distribution of area errors between H93 and GI old for all glaciers mapped in the same y of the glacier area showing Gaussian and student-t distributions fit to the data.The heavier tails distribution reflect the effect of larger blunders on glacier outlines, presumably due to varyi divides and debris cover/lateral moraine delineation between the area estimates.
Fig. 9. Distribution of area errors between H93 and GIold for all glaciers mapped in the same year, as a percent of the glacier area showing Gaussian and student-t distributions fit to the data.The heavier tails of the student-t distribution reflect the effect of larger blunders on glacier outlines, presumably due to varying hydrological divides and debris cover/lateral moraine delineation between the area estimates.
Fig. 9. Distribution of area errors between H93 and GI old for all glaciers mapped in the same year, as a percent of the glacier area showing Gaussian and Student t distributions fit to the data.The heavier tails of the Student t distribution reflect the effect of larger blunders on glacier outlines, presumably due to varying hydrological divides and debris cover/lateral moraine delineation between the area estimates.
timates by standard error propagation (root sum of squares) results in an error of 14 m a −1 .About 80 % of the observed length change rates in Epoch 1 and 2 (Fig. 8) are above this uncertainty.
All of the change parameters are sensitive to bias induced by interpretation uncertainty.In particular, the decision whether to include or exclude lateral moraines within the glacier area needs to be consistent within multi-temporal inventories.In Svalbard, glacier ice may exist beneath these lateral moraines (F.Navarro and A. Martín-Español, personal communication, 2013).Our inventories exclude lateral moraines by adopting a visual definition for delineating glaciers based on spectral appearance.Without widespread ground truth information about ice below debris, it is not possible to quantify potential error introduced.In addition, the decision whether to include or exclude lateral moraines is subjective and dependent upon the purpose of the glacier area outline.In Svalbard, the retreat of glaciers commonly occurs at the front rather than the sides; accordingly exclusion of the lateral moraines seems appropriate for studies of glacier extent changes.Alternately, using glacier area for volume estimation may require their inclusion (Radić and Hock, 2010;Huss and Farinotti, 2012;Martín-Español et al., 2013).

Discussion
Our new glacier inventory of Svalbard, GI 00s , can be used to extract spatial data reflecting topography and climatology of the archipelago.Median glacier elevation (Fig. 5a) is a characteristic of an individual glacier's hypsometry that is highly correlated with the equilibrium line altitude, or ELA (Braithwaite and Raper, 2009).The median elevation has been used for developing concepts of "glaciation limits" and proxies for the long-term ELA, the patterns of which suggest an inverse relation to the precipitation regime (Østrem, 1966;Schiefer et al., 2008).In Svalbard, spatial patterns of median glacier elevation have previously been used to infer spatial variability of the ELA (Hagen et al., 2003) and precipitation (Hisdal, 1985;Hagen et al., 1993;Winther et al., 1998;Sand et al., 2003).In central Spitsbergen, the spatial patterns of the median elevation match closely the mapped 1990 late summer snow line (a proxy for the ELA) distribution (Humlum, 2002).On Austfonna, snow depth variability shows a clear northwest-southeast gradient (Taurisano et al., 2007;Dunse et al., 2009) due to the predominance of precipitation coming from the Barents Sea (Schuler et al., 2007).This is reflected in lower median glacier elevations towards the southeast and higher towards the northwest (Fig. 5a).Spatial patterns of median glacier elevation similarly reflect the annual total number of melt days and summer melt onset as estimated from QuikSCAT scatterometry (Rotschky et al., 2011).The spatial patterns of median elevation over the archipelago reflect the local degree of glaciation, which is dependent upon both the terrain and the long-term regional climatological patterns that result in spatial variations of accumulation and ablation (mass balance) over the terrain surface.Higher glacier median elevations occur in the central drier regions of Spitsbergen and correspond to areas that experience lower average annual melt days, which could imply lower mass turnover.Moreover, these areas have lower percent glaciation and larger number of glaciers in the inventory (Fig. 4).
A glacier inventory represents a snapshot of the glacier geometrical extent, typically at a single point in time.Our glacier inventories are generally not a single point in time but cover a range of times, though each glacier outline has a distinct time stamp.Comparing multiple glacier inventories through time allows investigation of changes in some of the basic glacier geometry parameters.Changes in glacier area and length reflect the glacier's total response (Oerlemans, 2001).At smaller regional scales, area and length changes of individual glaciers manifest themselves differently to the presumably more or less uniform driving climate signal, due to variable glacier response times.Glacier response time is proportional to thickness and inversely proportional to the ablation rate at the terminus (e.g.Jóhannesson et al., 1989), such that front positions of small thin glaciers respond more quickly to the same climate change signal.Above a critical glacier size (i.e.larger glaciers) and holding all mass balance gradients similar, theory and modeling experiments predict a decreasing response time, with increasing glacier size resulting from the dynamic controls on response time (Bahr et al., 1998;Pfeffer et al., 1998).Estimated response times for Svalbard glaciers are in the range of decades to centuries, implying that observed front position changes still contain signals from earlier climatic events -especially true for Epoch 1. Epoch 2 changes may reflect climate changes during both Epoch 1 and Epoch 2. In Svalbard, the frequent surging behavior of many Svalbard glaciers (Lefauconnier and Hagen, 1991;Hagen et al., 1993;Hamilton and Dowdeswell, 1996;Jiskoot et al., 2000;Sund et al., 2009) complicates reconstructions of past climate from change records (as in Oerlemans, 2005).Finally, since each inventory represents a single glacier snapshot, the variation in temporal separation between inventories (epoch length), and its relation with the timing of each individual glaciers response will influence the observed average change rates.Enhanced interpretation of these changes may be possible with the inclusion of an accurate surge glacier inventory.This was not completed for this study due to insecurities in defining exactly which glaciers have fully surged or only partially surged (Sund et al., 2009) especially over the decadal time period of this study where glaciers that have surged may not be visible in the area changes given the long time span between the inventories.Future work should focus on generating such an inventory.
Figure 10 shows the difference between area/width retreat rates of Epoch 1 and Epoch 2 for the sample of 392 glaciers larger than 2 km 2 that have retreated in both epochs.About Average retreat rate di erence between epoch 1 and 2 (m a -1 ) Fig. 10.Annual average length changes calculated using the 'area/width' method for land-terminating (left) and tidewater (right) glaciers for Epoch 1 vs Epoch 2. Colors represent the difference in average retreat rates between the epochs.Black diamonds are glaciers that have advanced in either Epoch 1 or 2. Symbols in the maps (bottom) are identical to those in the scatter-plots (top).In the scatter-plots, points below the 1:1 line are glaciers that have experienced smaller retreat rates in Epoch 2 (blue) while points above the line are glaciers that have experienced greater retreat in Epoch 2 (red).Note the log scale of the scatter plot and linear color scale in the difference of average retreat rates.The size of the symbol is a log scale of the glacier size.34 Fig. 10.Annual average length changes calculated using the "area/width" method for land-terminating (left) and tidewater (right) glaciers for Epoch 1 vs. Epoch 2. Colors represent the difference in average retreat rates between the epochs.Black diamonds are glaciers that have advanced in either Epoch 1 or 2. Symbols in the maps (bottom) are identical to those in the scatter plots (top).In the scatter plots, points below the 1 : 1 line are glaciers that have experienced smaller retreat rates in Epoch 2 (blue), while points above the line are glaciers that have experienced greater retreat in Epoch 2 (red).Note the log scale of the scatter plot and linear color scale in the difference of average retreat rates.The size of the symbol is a log scale of the glacier size.
60 % of the differenced retreat rates are greater than the 95 % confidence interval of 10 m a −1 (see Sect. 4.4).For the entire sample and the significant subsample, ∼ 60 % of the glaciers have experienced larger retreat rates in Epoch 2, apparent as a shift in the histograms between the two epoch-averaged retreat rates (Fig. 8).The spatial patterns of Fig. 10 suggest a potential regional trend with greater retreat in Epoch 2 in southern Spitsbergen and less so in northwest Spitsbergen.Moreover, slight clustering or spatial autocorrelation in the retreat rate changes between the epochs is present.This may reflect the geometrical similarity between neighboring glaciers that is dictated by the topography and the more or less uniform driving climate signal.On top of the apparent clustering, larger outliers are present, with some neighboring glaciers having opposing (either positive or negative) differenced retreat rates.These local outliers represent variability in the individual glacier responses to the driving climate and/or the effect of past and present surge glaciers and their surge history in relation to the timing of the inventories.Combination with additional parameters such as geodetic volume changes (e.g.Nuth et al., 2007), estimated glacier volumes and thicknesses or terminus mass balance rates (Haeberli and Hoelzle, 1995;Hoelzle et al., 2007) will aid interpretations related to response times and climate.The most common parameter for inventory change in the literature is relative area change which is dependent on the glacier and terrain geometry (e.g.size, individual glacier tongue shape, bed topography etc.).Smaller glaciers often have larger relative changes, such that the variability with glacier size is heteroscedastic.This complicates statistical and spatial analysis of relative changes.To reduce the area dependency, we estimate length changes using 2 approaches (Sect.3.4).Area/width length changes are larger than the centerline length changes (Fig. 7a) as they integrate the entire front position change (including lateral changes), while centerline changes are one single measurement of the front.The area/width method minimizes area-related errors from uncertain upper glacier boundaries by limiting the analysis to area change below the median glacier elevation and may be fully automated allowing easy retrieval from repeat inventories.
Using the entire sample of glaciers that exist consistently in GI old , GI 90 and GI 00s , the total area change rates and relative area change rates are 22 % smaller in Epoch 2 (−24 km 2 a −1 , 0.23 % a −1 ) than in Epoch 1 (−31 km 2 a −1 , 0.26 % a −1 ).Alternately, summing the area/width (centerline) length changes for all glaciers results in 14 (30) % more negative length changes during Epoch 2 (−29.5 (−15.3)km a −1 ) as compared to Epoch 1 (−25.6 (−11.9)km a −1 ).Thus, while area change rates were larger in Epoch 1, the length change rates were larger in Epoch 2, which is similarly reflected in Fig. 7b.The total summed change rates of glacier tongue width during Epoch 1 is −2.2 km a −1 as opposed to −1.4 km a −1 for Epoch 2. General reduction in glacier area loss rates in Epoch 2 is probably a geometrical effect of decreased lateral wastage of the glacier tongues, despite the faster average retreat rates experienced in Epoch 2.

Conclusions
This study describes the creation of a consistent multitemporal digital glacier inventory of the Svalbard archipelago, based on the structure of the previous inventory (Hagen et al., 1993).Our new digital inventory is based on historic data that are available digitally and then progressively updated through time to maintain consistency between the glacier outlines.This required modification of the identification system already in place (WGMS, 1989;Hagen et al., 1993) for glaciers that have retreated and separated.Moreover, the newest inventory also includes snowpatches and glacierets that are less than 1 km 2 as identified in available cloud-free SPOT, ASTER and Landsat images.GI 00s , the present digital inventory, coheres to both GLIMS and WGI standards (with slight modifications) and has been incorporated into the Randolph Glacier Inventory (Arendt et al., 2012) and the GLIMS database.The inventories may also be downloaded from the Norwegian Polar Institute data archive (http://data.npolar.no/dataset/89f430f8-862f-11e2-8036-005056ad0004).
In total, the GI 00s inventory of glaciers in Svalbard contains 1668 individual glacier units, for a total glacierized area of 33 775 km 2 , or ∼ 57 % of the archipelago.About 60 % of the glacierized area is on Spitsbergen and 40 % on the island to the east-northeast.Between 168 and 197 tidewater glaciers (depending on how a single tidewater glacier is defined) drain 68 % of the glacierized area through a summed glacier terminus width of about 740 km.The glacierized area has decreased by an average of ∼ 80 km 2 a −1 over the past ∼ 30 yr, a reduction of 7 %.For a sample of ∼ 400 glaciers in south and west Spitsbergen, glacier retreat was greater after 1990 (Epoch 2), while area change rates were greater in the decades before 1990 (Epoch 1), corresponding to more lateral wastage in the early period.
We suggest that reducing the dimensions of area change to a length scale may provide a more useful parameter for spatio-temporal analysis of change signals.This falls in line with previous investigations (Raup et al., 2009) and may enhance the use and incorporation of glacier inventory changes, for example, into numerical models (e.g.Oerlemans, 1997;Vieli et al., 2001;Nick et al., 2009) and/or temperature reconstructions (Oerlemans, 2005;Leclercq and Oerlemans, 2012).The spatio-temporal variability of the length change rates suggests response time variation that requires further investigation.With increased accuracy and capability to automatically generate glacier inventories at higher temporal resolutions using satellite data, the generation of repeat glacier inventory changes in the form of area, tongue width and length changes will become an important observational dataset for future glacier-climate studies.frequency of the bumpy glacier surface seen in Fig. A1a.A single frequency of the artifacts was not practically determinable for use in a frequency stop filter due to the similarity of the artifacts with the natural glacier topographic fluctuations and the frequency of off-glacier terrain.Therefore, we choose to apply a standard low-pass filter in the frequency domain, using a Hanning window.The set of parameters (order and frequency cut-off) is chosen by minimizing the standard deviation of the glacier differences between the GDEM and ICESat.Post-filtered topography was not sensitive to small fluctuations in the parameters.The filtered GDEM (Fig. A1b) shows improvement over smooth glaciers (Fig. A1e) but resembles a downgraded (lower resolution) product over the rougher surrounding terrain (Fig. A1f).Most of the higher-frequency noise is removed (Fig. A1g), though the lower-frequency bumps remain, with maximum differences of up to 50-60 m (Fig. A1h).Therefore, our final post-processed GDEM is a compilation of the Fourierfiltered glacier surface with a median block filtered nonglacier surface.

Fig. 1 .
Fig. 1.Spatial and temporal coverage of glacier inventories (GI).The three maps (top) and the filled bars (bottom) show digitally available outlines.The black unfilled bars are the H93 tabular inventory (Hagen et al., 1993).The satellite-based inventory, GI 00s (in red), is the first digital inventory covering the entire archipelago.The shades of individual colors in the maps represent time within each inventory.Lighter shades are earlier; i.e. the lightest green is from 1936 and the lightest red from 2001.Grey represents no data.

Fig. 2 .
Fig.2.Regions, primary and secondary drainage basins (i.e. the first three digits of the IDs).The colors depict the primary drainage basins (i.e. the first two digits of the IDs) and correspond to the colors in the glacier hypsometries shown to the right.Glacier hypsometries are extracted from the ASTER GDEM.

Fig. 3 .
Fig. 3. Examples of the digital inventories from a selection of (a) land-terminating and (b) tidewater glaciers in southern Spitsbergen.Selected glaciers are shown with their local identification codes, centerlines and glacier tongue widths.Note that only centerlines and tongue widths from the lowest 10% of the centerline of the GI 00s are shown.

Fig. 4 .
Fig. 4. (a) Glacier number and area distribution of the GI 00s inventory.E.g. 80 % of the glaciers in the inventory are less than 10 km 2 , which makes up less than 10 % of the total Svalbard glacierized area.(b) Percent glacierized area of the secondary drainage basins.(c) Number of glaciers within secondary drainage basins reflects an inverse relationship to the percent glacierized area.(d) Percent tidewater glacier area for each secondary drainage basin shows the dominance of tidewater glaciers in Nordaustlandet and the three ice field clusters in Spitsbergen (south, northwest and northeast).

Fig. 5 .
Fig. 5. Spatial distribution of (a) median elevation, (b) mean slope in degrees and (c) mean aspect in degrees from north.The color scale is shown as a histogram with the number of glaciers.The distribution of mean slope is slightly bimodal, reflecting two styles of glaciation: the flatter glaciers and icecaps that fill valley floors, and the steeper cirque-style glaciers that sit higher up along the valley walls.The distribution of mean aspect per glacier suggests a dominance of north-facing glaciers; however, histograms of all glacier DEM pixel aspects show a uniform distribution with aspect.The histogram thus reflects the dominance of small glaciers facing northward.

Fig. 6 .
Fig. 6.(a) Distribution of land-terminating and tidewater glaciers in Svalbard.The sum of tidewater glacier front widths within each secondary drainage basin is shown as size (area)-proportional circles.(b) Relative glacier area changes between H93 and GI 00s for each secondary drainage basin, with total glacier area shown as size (area) proportional symbols.

Fig. 7 .Fig. 7 .
Fig. 7. (a) Temporal mean retreat rates for Epoch 1 and 2, as estimated from the centerline differences and by the area change divided by tongue width ("Area/width").The inset shows the histogram of differences between the two retreat rate estimates.The "Area/width" retreat rates are larger than the centerline due to both incorporation of lateral losses in the area-width retreat rate estimate and a lack of representative retreat sampling of the centerline.(b) Relative changes during Epoch 1 and Epoch 2 for area changes, centerline length changes and the "Area/width" length changes.

Fig. 8 .Fig. 8 .
Fig. 8. Average length change rates calculated using the "area/width" method between (a) GIold and GI90 and (b) GI90 and GI00s.The inset map is a section of northwest Spitsbergen.The histogram insets show the number of glaciers for each color used in the map.The straight bar line in the histograms are of the alternate epoch for comparison.The distribution of Epoch 2 length changes is more negative than that of Epoch 1. 28

Fig
Fig. A1.(a) Hillshades draped by elevation for the original ASTER GDEM, (b) the Fourier-filtered GDEM, and (c) a 2008 SPOT5-HRS DEM.(d) The differences between the filtered GDEM and the SPOT DEM with the black areas containing correlations that are less than 40.(e) and (f) show normalized histograms (y axis as percentage of pixels) of elevation differences (x axis in meters) for land and glaciers, respectively, for the DEMs and ICESat.In the legend of (e), "GDEMo" and "GDEMf" are the original and filtered GDEM, respectively.(g) and (h) show two individual ICESat profiles over higher-elevation glacier regions (typical for low visible contrast areas) from 27 May 2005 and 19 March 2008, respectively, with their differences to the three DEMs.The thick brown line shows the topography, while the other colored lines correspond to the labels presented in (e).

Table 3 .
List and description of attributes available for each glacier inventory.