Glacier area and length changes in Norway from repeat inventories

Introduction Conclusions References


Introduction
Glaciers are key indicators of climate change, making their monitoring important (e.g.Vaughan et al., 2013).Remote sensing techniques are ideal for measuring glaciers on a large scale, as they cover remote glacierized areas with relatively little effort.Optical images provided by the Landsat TM/ETM+, Terra ASTER, or SPOT HRS have proven to be very efficient for mapping glacier extents (e.g.Paul et al., 2002;Paul and Kääb, 2005;Paul et al., 2007;Bolch et al., 2010;Nuth et al., 2013).
Glacier outlines are typically obtained from satellite images using thresholded multispectral band ratios (Bayr et al., 1994;Sidjak and Wheate, 1999;Paul and Kääb, 2005;Kargel et al., 2014).An important advantage of the Landsat TM/ETM+ sensors is their large swath width, which is ideal for mapping extensive glacier regions.
Glacier inventory data are used for modelling glacier mass balance (e.g.Marzeion et al., 2012;Radić and Hock, 2013), estimating ice volumes (e.g.Huss and Farinotti, 2012;Grinsted, 2013;Andreassen et al., 2014), or predicting global sea level rise (e.g.Leclercq et al., 2011;Gregory et al., 2013).Despite Norway's long tradition of monitoring glaciers, there are still few data available on spatiotemporal change.Due to Norway's favorable topography and climate, hydro-power accounts for 98 % of its electricity; about 15 % of the run-off comes from watersheds that are partly glacierized (Andreassen et al., 2005).For this reason, the Norwegian Water Resources and Energy Directorate (NVE) maintains a database of glaciological data in mainland Norway.Nevertheless, most of the glaciers have not been mapped due to their remote locations.In addition, most glaciers in Norway lacks information about the spatial and temporal variations of glacier change.
Previous glacier inventories of Norway from the 1960s to the 1980s lack digital glacier outlines (Østrem and Ziegler, 1969;Østrem et al., 1973, 1988), and glacier area and length assessments are complicated e.g., by unknown ice divides.The most recent satellite-derived glacier inventory of Norway is based on Landsat TM/ETM+ (Andreassen et al., 2012b).It uses a GIS-based approach and is compiled following the Global Land Ice Measurements from Space (GLIMS) guidelines (Kargel et al., 2005;Racoviteanu et al., 2010).This data set is a highly detailed digital baseline product ideal for glacier area and length change assessments (Andreassen et al., 2008;Paul and Andreassen, 2009;Paul et al., 2011).Glacier length measurements are considered one of the most important ways to quantify glacier change in the future (Hoelzle et al., 2003).Historical length change observations can give much information about how glaciers respond to climate (Leclercq et al., 2014).The newly compiled Norwegian glacier inventory is available through the GLIMS database and as a published book (Andreassen et al., 2012b).Globally, areas with multi-temporal glacier outline data sets are rare (Kargel et al., 2014).Because long term glacier change assessments are crucial for understanding glacier response to climate (Hoelzle et al., 2003), there is a need to complete such a multi-temporal glacier inventory.
For the first time, we present multi-temporal data sets derived from Landsat TM/ETM+ satellite images and topographic maps for all glacierized areas in Norway.We perform a glacier area and length change assessment which is based on three data sets from the time ranges of 1947 to 1985 (GI n50 ), 1988 to 1997 (GI 1990), and 1999to 2006(GI 2000 ).We compare in situ length change observations with length changes from automatically derived centerlines.We extend the data sets prior to GI n50 using older topographic maps, which allows us to conduct an extended glacier area and length assessment on five ice caps in northern Norway.Concerning the multi-spectral band ratio technique, we demonstrate that mapped glacier areas are sensitive to small variations in the chosen ratio thresholds.

Study region
Mainland Norway extends from 58 • to 71 • N and from 5 • to 31 • E and covers an area of 385 199 km 2 (Fig. 1a).The identified 2534 glaciers in the most recent glacier inventory have a total area of 2692±81 km 2 (using ±3 % as uncertainty), covering 0.7 % of the area of Norway (Andreassen et al., 2012b).In the most recent glacier inventory, glacier complexes are divided into individual glacier units.These glacier units share common divides if they are part of a glacier complex, otherwise they correspond to single glaciers without a drainage divide.The number of glacier units in the most recent glacier inventory is 3143.We divided the study area into three geographical parts: northern, central, and southern Norway, which were further split into 36 sub-regions (map in Figs.2e).
Coastal regions in Norway have a warm and moist maritime climate, while the interior is drier and colder.Climate gradients along a west-east transect are pronounced, especially in southern Norway.This west-east pattern is caused by the westerly winds and the Gulf Stream, together with the shading effect in the eastern parts due to the coastal mountains (Hanssen-Bauer et al., 2009).These climatic factors contribute to warmer conditions in Norway compared to similar latitudes elsewhere in the world.Norway has a latitudinal gradient in terms of mean temperature and precipitation, which both decrease from south to north.However, along the coast, there is no pronounced variation on climate because of the ice-free Norwegian Sea, although Norwegian glaciers span over ∼ 1500 km from south to north (Fig. 1b).The mean equilibrium-line altitude (ELA) of the glaciers increases inland, and decreases towards the north due to climatic differences (Andreassen et al., 2005).
The first systematic glacier length observations in Norway were initiated around 1900 (Hoel and Werenskiold, 1962).Throughout the 20th century, Norwegian glaciers have generally retreated, except for intermittent advances of the coastal glaciers.Periods of increased winter precipitation have contributed to a temporary mass gain for all glaciers.Advances were recorded around the years 1910 and 1930, and in the 1990s (Andreassen et al., 2005;Nesje et al., 2008).Since the beginning of the 2000s, all glaciers monitored by NVE have been in a state of retreat (Andreassen et al., 2005;Winkler et al., 2009).
Glacier inventories of Norway were published in 1973 for northern Norway (Østrem et al., 1973) and 1969 and 1988 for southern Norway (Østrem and Ziegler, 1969;Østrem et al., 1988).The first complete and satellite remote sensing-based inventory of Norway was published in 2012 (Andreassen et al., 2012b).In this paper, Norway refers to mainland Norway only.Area and length changes for Svalbard were recently published by Nuth et al. (2013).
3 Data and methods

Data set background
Our glacier inventory data are compiled using multi-spectral Landsat satellite data for the time periods GI 2000 and GI 1990 , topographic maps based on aerial photographs for GI n50 , and analogue maps to extend glacier outlines further back in time (GI 1900 , prior to the GI n50 data set) (Fig. 2).In our analysis, we compare the data sets resulting in three epochs: full epoch (GI n50 -GI 2000 ), epoch 1 (GI n50 -GI 1990 ), and epoch 2 (GI 1990-GI 2000 ).
In epoch 1 and epoch 2, some glaciers had less than 10 years between the two data sets compared, corresponding to 12 % of the numbers of glaciers in both GI 1990 and GI 2000 (Table 1).Ideally, glacier inventories should be retrieved in intervals of a few decades when used in change assessments, to account for the glacier response time (Haeberli, 2004).However, if a glacier region encounters very fast down-wasting of the glaciers, shorter mapping intervals can be used, which is the case for many Norwegian glaciers.
The multi-temporal data sets contribute to the monitoring of essential climate variables (ECVs) and follow the Global Climate Observing System principles (GCOS, 2003).The data sets were created in accordance with the guidelines on how to monitor glaciers and ice caps established by the Global Terrestrial Network for Glaciers (GTN-G, 2009).

GI 2000 and GI 1990 -Landsat satellite imagery
The Landsat TM/ETM+ satellite images have multiple advantages compared to imagery from ASTER and SPOT due to: 1) the larger swath width of Landsat, 2) better availability of Landsat images, as other optical satellites were not operational during the time periods, and 3) Landsat has freely available georeferenced and orthorectified satellite scenes.The year of satellite acquisitions and the spatial coverage for the GI 1990 and GI 2000 data sets are presented in figure 2c and d.GI 1990 and GI 2000 span over periods of 9 and 7 years respectively, as it proved impossible to map outlines for all Norwegian glaciers within one or a few years using Landsat TM/ETM+.This is due to a lack of cloud-free Landsat TM/ETM+ satellite scenes as a result of Norway's pronounced maritime climate.Seasonal snow cover, due to the high precipitation rates throughout all seasons, also makes satellite image interpretation challenging (Andreassen et al., 2008).Due to extensive cloud coverage and partly also seasonal snow, full coverage for the GI 1990 was not possible.No usable scenes were available for Jostedalsbreen, Lofoten/Hamarøy and part of inner Troms (see Sect. 3.8).
Prior to the derivation of glacier outlines from the Landsat scenes, we carried out an accurate orthorectification and quality check of the images using PCI Geomatica © .The Landsat L1T/L1Gproducts were delivered orthorectified and often used as-is after a quality check (Table 2).However, selected satellite images had to be orthorectified prior to the derivation of outlines due to insufficient quality of the L1T/L1G-products, especially in mountainous areas.The root mean square values (RMSE) for both the purchased satellite scenes and the orthorectified ones had an accuracy of less than ∼ 30 m.We calculated the band ratios for the Landsat images by including the red band (TM 3 ) and the short wave infrared band (TM 5 ).We decided to use TM3 TM5 instead of TM4 TM5 for the glacier delineation following Andreassen et al. (2008) results from glacier outlines in Jotunheimen.They show that TM3 TM5 performed better for ice located in shadow and for debris-covered ice compared to TM4 TM5 .The band ratio method uses threshold values optimized for each satellite scene.We used TM3 TM5 ≥ t 1 , where t 1 varied between 1.6 to 2.8.To improve results in shadowed areas, we included an additional threshold on the blue band, TM 1 ≥ t 2 , where t 2 is either 35 or 60, with some exceptions (Table 2) (Paul and Kääb, 2005;Paul and Andreassen, 2009).We applied a median filter on the glacier outlines to eliminate individual glacier pixels.Outlines were further manually corrected in case of debris cover, glacier-lake interfaces, clouds or cast shadow which hampered the automatic mapping.Only very few outlines had to be corrected for debris cover since the glaciers in Norway are mostly debris-free.Lakes and seasonal snow misclassified as glaciers were masked out from the glacier outline product.We used topology editing in ArcGIS © for the manual corrections and the delineations of ice divides.Topology rules allowed features that share the same geometry to be updated simultaneously.The methods of filtering, human inspection and editing of the data sets are described in the glacier inventory by Andreassen et al. (2012b).

Band ratio accuracy and threshold sensitivity
The accuracy of the band ratio method and the sensitivity of the used threshold values are essential for change assessment of glaciers.Orthophotos from the same acquisition year as the satellite images are ideal for determining accuracy, but are rarely available.In Jotunheimen, a mountainous region in southern Norway, glacier outlines were compared with orthophotos taken one year apart, and an area difference of −2.4 % was found (Andreassen et al., 2008).Fischer et al. (2014) show that Landsat derived outlines (year 2003; medium spatial resolution (30 m)) compared to orthophotos (year 2003; high spatial resolution (50 cm)) for eastern Switzerland show similar results meaning there is comparable accuracy between the medium-resolution and high-resolution source data for glaciers > 1 km 2 .On the other hand, they found that glaciers < 1 km 2 , the uncertainty of the outlines increased with decreasing glacier size.For debris-free glaciers, the band ratio method is robust and accurate (Albert, 2002;Paul et al., 2003) with an accuracy between ±2-5 % (e.g.Paul et al., 2013).
Here, we operate with an accuracy of ±3 %, implying that the inventory of Norwegian glaciers has a total accuracy of 2692±81 km 2 (Andreassen et al., 2012b) (2668±80 km 2 for the glaciers included in GI 2000 ).The automatic band ratio method and manual digitizations give similar results, but the band ratio method often obtains a smaller glacier area as it tends to exclude some mixed pixels (Paul et al., 2013).
The mapped glacier area depends strongly on the chosen threshold value.The sensitivity of selected threshold values used on the ratio TM3 TM5 ≥ t 1 and the additional blue band TM 1 ≥ t 2 have been investigated for 57 glacier units in western Finnmark, Northern Norway.A Landsat 5 TM satellite scene with good snow and cloud conditions from the year 2006 was used (Area code 1 2000 in Table 2).By calculating the difference in number of pixels mapped for selected thresholds, a percentage difference of area relative to the applied threshold is calculated (Fig. 3).We used TM3 TM5 ≥ 2.4 and TM 1 ≥ 35 or 60 as reference threshold values (yielding a total of 69.3 and 65.6 km 2 respectively).

The TM3
TM5 ratio thresholds range from 2.0 to 2.8 with increments of 0.2.There were some outliers strongly affecting the mean values of area change between the thresholds compared, thus it is more representative to use median values (Fig. 3a).
Comparing the area derived from the thresholds TM3 TM5 ≥ 2.0 and 2.4, and TM 1 ≥ 35 yielded a median area increase of 12 %.With this change in threshold a larger glacier area was mapped compared to the reference threshold value (also glaciers in cast shadow).With lower threshold values more noise were included i.e. mixed pixels containing snow/ice and rock/debris.Similarly, when comparing TM3 TM5 ≥ 2.4 and 2.8, and TM 1 ≥ 35 we find a median decrease in area of −11 % (−3.1 km 2 ).The higher threshold values used for TM3 TM5 reduce noise, but include less glacier area compared to lower threshold values, due to less mixed pixels including both ice and terrain features.The TM3 TM5 threshold should be as low as possible to include the dirty ice around glacier perimeters (Paul et al., 2013).

If TM3
TM5 ≥ 2.4 was used with TM 1 ≥ 60, we found less variation when varying the threshold values compared to using the TM 1 ≥ 35.This means a median area decrease of −4 % (−1.2 km 2 ) using TM3 TM5 ≥ 2.4 to 2.8 and a median area increase of 3 % using TM3 TM5 ≥ 2.0 to 2.4.By applying a threshold of TM 1 ≥ 35, the area mapped was more sensitive and included more mixed pixels compared to TM 1 ≥ 60 (Fig. 3b and c).Glaciers in cast shadow were most sensitive to variations in the applied threshold and manual onscreen digitization was necessary in some occasions.

GI n50 -topographic maps
The GI n50 data set was derived by digitizing 168 first edition 1 : 50 000 topographical maps in the N50-series from the Norwegian Mapping Authority based on aerial photographs acquired between 1947 and 1985 (Fig. 2b).Digital form is required to conduct change analysis and to compare glacier drainage basins using the same drainage divides.(Andreassen et al., 2008).
The first edition N50-paper maps were scanned and georeferenced using ground control points in a reference map (from European Datum zone 32, 33, and 34 to WGS 84 UTM zone 33).The glacier outlines were then digitized on-screen.We used a first-order polynomial transformation and obtained RMSE values of less than 10 m.The years of the aerial photographs were provided by the Norwegian Mapping Authority, and the acquisition years were also checked on every map sheet.Some map sheets (e.g.1532-2 Altevatnet) have several and thus unknown exact acquisition years (e.g. between 1947 and 1951).In those cases, the first mapping year (e.g.1947) was allocated to the glaciers within the map sheet.
Three investigators independently digitized a sample of 10 glaciers in the Frostisen region ranging from 0.06 to 0.93 km 2 by size (Frostisen map sheet 1331-2) to estimate the uncertainty arising from the subjective interpretation by the digitizer.65 % of all glaciers in the GI n50 data set are within this size class.The standard deviations of the total glacier area were in the range of 0.0018 to 0.0066 km 2 for the 10 selected glaciers digitized by three different investigators.This indicates little variation in the digitizing accuracy between the interpreters and minor effects on the end result.We expect a larger uncertainty for the derivation of the actual outlines on the maps.However, an uncertainty analysis was not applied to the topographic maps, because the orthorectification and digitizion of the imagery would be a huge amount of work.Another uncertainty factor was the unknown working methods and mapping principles of the cartographers of the time.

GI 1900 -analogue maps
The applicability of older analogue maps was tested in the two northernmost glacier regions in Norway (Seiland and Øksfjord, see Fig. 2a).A map series named Gradteigskart was used, which are quadrangle maps constructed from land surveys (map scale 1 : 100 000).The surveys of Seiland and Øksfjord were conducted in the period 1895-1907(GI 1900 ) and cover three map sheets (174-92:1895, 174-106:1896, and 147-178:1907).The map sheets include the five largest ice caps in the two regions, Nordmannsjøkelen, Seilandsjøkelen, Øksfjordjøkelen, Svartfjelljøkelen, and Langfjordj økelen.The GI 1900 -map sheets were scanned and georeferenced using local transformation methods.The glacier outlines were manually digitized from the georeferenced maps.Glacier outlines from the GI 1900 data set are less accurate than GI n50 .Old analogue maps can have severe planimetric distortions due to complex topography, and in certain cases, glacier extents are known to be overestimated (Østrem and Haakensen, 1993).While it is worthwhile to incorporate these outlines into the change analysis, the results must be interpreted with care and considered as an estimate rather than an accurate benchmark.
For three composite glaciers in West-Finnmark (Langfjordjøkelen, Øksfjordjøkelen, and Svartfjelljøkelen), we tested four transformation methods for georeferencing: 1) spline, 2) adjust, 3) second-order polynomial, and 4) third-order polynomial.We used 25 ground control points in each map sheet with a mean total RMSE of ±50 m.We chose to use the second-order polynomial as transformation method, as it had the best agreement with real topography and the other data sets.In total, the area of the second-order polynomial transformation method was 89 km 2 .Area differences between the tested methods varied up to 1.8 % (1.6 km 2 ) relative to the applied method.

Digital Elevation Model (DEM)
A DEM is required for deriving topographic parameters, ice divides, and centerlines.We used GIS analysis to extract glacier parameters from the data sets.Topographic parameters were calculated for GI 2000 , including minimum, maximum, mean, and median elevation, aspect, and slope.The spatial resolution of the DEM provided by the Norwegian Mapping Authority is 20 m, and the elevation contours, which the DEM was constructed from, were derived from aerial photography.The DEM from the Norwegian Mapping Authority was used for all ice divides, except on the Folgefonna and Hardangerjøkulen ice caps where DEMs based on LIDAR data acquired in 2007 and 2010 were available (Andreassen et al., 2012b).The national DEM was used to calculate the topographic parameters for all glaciers.Accordingly, the date of the DEM is not always coincident with the date of the outlines.In the case of GI 2000 , the DEM data are up to ∼ 20 years older than the outlines.Still, unless the elevation changes are very large, the surface changes have only a minor impact on many of the derived parameters due to the extensive number of pixels averaged when calculating slope, mean altitude, and aspect (Frey and Paul, 2012).However, minimum altitude is sensitive to the applied DEM, as the minimum elevation is highly affected by rapid changes in the glacier termini.Ideally, we would have had multi-temporal DEMs, with each set of outlines having its own DEM, but this is not available for our study area.Paul et al. (2011) explored the possibility of using the ASTER global DEM in the Jostedalsbreen glacier region.The study revealed too many artifacts close to the termini of the glaciers, and concluded that the national DEM was a better choice.

Division of glaciers
To create an inventory of individual units, glaciers were divided into glacier units based on ice divides.We used the Hydrology tools Flow accumulation and Flow direction in ArcGIS © .Flow accumulation creates a raster of accumulated flow into each cell, and Flow direction finds the steepest down slope neighbor for each cell.Both data sets are created from the national DEM.Additionally, discharge basins from the NVE were used to define ice divides manually.In the GI 1900 , GI n50 and GI 1990 , each glacier complex was separated using the same ice divide as GI 2000 .The glacier basins were adjusted to include glacier outlines from all three data sets.If GI 1900 , GI n50 , and GI 1990 extended the GI 2000 perimeter, e.g, due to disintegrating glaciers, the ice divides and glacier basins were adjusted using topology editing and the DEM.In GI 2000 the 36 glacier regions were arranged from north to south (map in Figs.2e), and within the regions, glacier IDs (1-3143) were assigned automatically from north to south based on latitude.The Norwegian Glacier Inventory book includes maps and tables for all mapped glacier entities covered in GI 2000 (Andreassen et al., 2012b).

Deriving centerlines
Glacier length is reported for all glacier units in the previous Norwegian glacier inventories (Østrem and Ziegler, 1969;Østrem et al., 1973, 1988).However, in the latest inventory of Norwegian glaciers, glacier length was not included (Andreassen et al., 2012b).In this study, we extend the inventory by adding glacier length.
We calculated centerlines by applying a three-step cost grid -least cost route approach which requires glacier outlines and a DEM as input (Kienholz et al., 2014).We excluded glaciers smaller than < 1 km 2 (from GI n50 ) to reduce the noise from seasonal snow cover, which gives relatively large errors over smaller areas.Our method for deriving centerlines was reproducible and fast compared to digitizing flow lines by hand, which would be very time consuming.In the first step, the algorithm identifies the lowest glacier points as termini, while using local elevation maxima along the glacier outlines to identify heads for each major glacier branch.In the second step, the algorithm determines the actual centerlines between heads and termini by calculating the least cost route on a cost grid (Fig. 4).The cost grid is established using the grid cells' elevations and Euclidean distances to the glacier edge.To allow for plausible centerlines, the lowest cost values are allocated to the glacier center and to lower glacier reaches.In the third step, the algorithm divides the centerlines into individual branches, which are then classified according to a branch order (Kienholz et al., 2014).Each step of the algorithm is implemented independently so that results can be checked visually and corrected if necessary.
For each glacier, the longest centerlines were extracted to describe glacier length, in agreement with previous studies (e.g.Paul et al., 2009).To calculate length change, we took the longest centerlines from the GI n50 dataset as a reference and clipped them with the more recent glacier outlines in the frontal areas (Fig. 4).We visually checked whether the clipped centerlines ran through glacier termini in GI 1900 , GI 1990 , and GI 2000 .Edits were required in case of glacier advance, or perpendicular shift of the termini relative to the original centerline.In these cases, the centerlines were manually lengthened following contour lines.Manual corrections were particularly important for the dataset GI 1900 due to larger glacier areas than in GI n50 .The length change was derived from the length differences of the centerlines.Tests indicated that this procedure yielded more plausible length changes than recalculating the centerlines for each period using the corresponding outlines.
The main problem with such recalculation was that glacier recession often leads to the emergence of nunataks.The centerlines run around these new nunataks, increasing the glacier lengths, which interferes with the goal of isolating the change of the termini positions only.
In many cases the comparison between the centerlines from the glacier inventories was not feasible.For example because GI 1900 did not always align with the other data sets, preventing measurements of the length change.In case of Nordmannsjøkelen, only 1 out of 9 glacier units could be used (Table 3).
Cumulative time series of glacier front position measurements were available from the NVE database, and we compared these in situ length change measurements with our length changes derived from centerlines for the full epoch.We used data from 12 glaciers with corresponding measuring periods in both the in situ measurements by NVE and the cumulative length changes from topographic maps and satellite images derived in this study (Table 4).Additionally, five glaciers were compared for epoch 1, and six glaciers for epoch 2. Some of the in situ measurements began before or after the GI n50 first mapping year, but series were included if the gap was no larger than 5 years.
Using an uncertainty of ±1 pixel for satellite images, variations less than 30 m can not be identified (Paul et al., 2011).

Inclusion, exclusion, and representation of data
In order to provide a complete picture of glacier variations in Norway, the objective was to include as much data as possible in the analysis.However, due to the heterogeneous nature of the data sets, parts of the data for different analyses were excluded .From the GI n50 and GI 1990 data sets, only glaciers that overlap in space with the GI 2000 data set were included.This approach was chosen because of the two earlier data sets tend to have higher uncertainties.For our analysis, 400 snow-ice patches was included that could be remnants of glaciers in the GI 2000 glacier areas, to make a more precise analysis of the area change.We assumed the snow fields were remnants of glaciers if they were located within previous glacier outlines older than GI 2000 .This results in a total of 2722 glacier units, an area of 2668 km 2 in the glacier area change (GAC) analysis for the full epoch divided into 1396 glacier units in southern Norway, 666 glacier units in central Norway, and 660 glacier units in northern Norway.For epochs 1 and 2, 1684 and 1953 glacier units were included in the analysis, respectively.In the glacier length analysis, we included 564 glaciers for the full epoch, 286 for epoch 1, and 283 for epoch 2. The reason for the significantly lower number of glaciers in the analyses for epoch 1 and 2 compared with the full epoch is insufficient satellite imagery due to cloud cover and seasonal snow for the dataset GI 1990 (see Sect. 3.2).Fig. 2c presents the areas where data were lacking for GI 1990 .For example, during the Landsat TM 4/5 and ETM+ 7 operation period , the largest ice cap, Jostedalsbreen, in western Norway where most of the data are missing, had only one Landsat satellite acquisition with preferable mapping conditions (from 16 September 2006).
The mean area change for glaciers not overlapping in all epochs was −0.119 km 2 (2722 glacier units).Including only glaciers overlapping in all epochs, mean area change was −0.096 km 2 (1684 glacier units).This variation indicates that Jostedalsbreen and glaciers in western Norway have been retreating.
In many cases, change assessments were challenging due to adverse snow conditions in the GI n50 data set.Visual inspection was done to identify snow patches often accumulated in gullies and ridges.
We used glacier basins to cut perennial and seasonal snow attached to the glaciers.In total 251 glacier units from the GI n50 data set were split using the glacier basins, and parts that were assumed to be seasonal snow (119 km 2 ) were detached and removed.All in all, 396 km 2 was excluded from the GI n50 data set, including parts of glacier units and additional single snow patches.In GI 1990 , 64 km 2 was excluded outside the basins, and 183 glacier units were cut using the basins.However, the numbers are not comparable due to incomplete coverage of glaciers in the GI 1990 data set.

Glacier area changes
The glacierized area in mainland Norway has decreased from 2994 km 2 in GI n50 to 2668 km 2 in GI 2000 , an area reduction of −326 km 2 (−11 %) (Table 5).Although most glaciers were subject to area loss, 20 % of glaciers increased in area during the full epoch (in total 37 km 2 ).For epoch 1, 19 % of glaciers increased in area, corresponding to 7.4 km 2 .In epoch 2, 63 % of glaciers gained area, corresponding to 98 km 2 .All in all, we find a clear decrease in glacier area for the full epoch (−10.9 %) as well as for epoch 1 (−10.5 %), but an increase in glacier area for epoch 2 (2 %) (Table 5).However, epoch 2 shows an area gain for the northern (4.7 %) and central parts (8.3 %), whereas the southern region decreased in area (−7.1 %).
In Fig. 5, we present glacier area change for northern, central, and southern Norway for the full epoch using normalized values ( GAC root(A) where A is the initial glacier area (Raup et al., 2009).This allows us to compare different groups without exaggerating the influence of the small glaciers, as is the case with values given in percentage, or the large glaciers, in the case of area changes given in absolute values (km 2 ) (Raup et al., 2009) (see Sect. 4.5).In Figs.5b and c, we find less glacier area decrease in the central part of Norway (glacier regions 13-19) for the full epoch.In the northernmost glacier regions 1-4, we find the strongest retreat pattern of the Norwegian glaciers (Fig. 5b).This is in line with in situ observations from the only NVE-monitored ice cap in this area, Langfjordjøkelen (Glacier region 2), which shows a strongly negative mass balance and area reduction over the last few decades (Andreassen et al., 2012a).The trend of retreating glaciers has also been seen on Svalbard, north of mainland Norway.Svalbard's glaciers show a total area reduction of −7 % the last ∼ 30 years (Nuth et al., 2013).
Glaciers in northern Norway are located at lower elevations than glaciers in southern Norway, which may explain their stronger retreat seen over epoch 1 and the full epoch (Fig. 5 and Table 5).
The distribution of area with elevation is presented in Fig. 6 for northern and southern Norway for the GI n50 and GI 2000 data sets, illustrating glacier area changes with elevation for the full epoch.
Northern and central Norway are grouped because they show similar area-elevation distributions.
For both southern and northern Norway, area decrease was larger at lower elevations than of high elevations, with area changes observed at all elevations (Fig. 6).In the elevation range between 1000 to 1700 m a.s.l., corresponding to 62 % of the total elevation range, the total absolute area loss is 201 km 2 .Many of the largest ice caps in both southern and northern Norway are located in this elevation range, and there are only a few glaciers above 1700 m a.s.l. in central and northern Norway.

Glacier length changes
The total centerline length for all glaciers in the GI 2000 data set (including all 3143 glacier units) is 3282 km thus, average glacier length is 1 km.The average length change for the full epoch is −240 m (Table 6).For comparison, we chose 286 glaciers that were included in all the data sets and whose area is > 1 km 2 in GI n50 .In this group of glaciers, the mean glacier centerline lengths are very similar between the two recent data sets, 2.86 km for GI 2000 and 2.91 km for GI 1990 , while the GI n50 has a longer mean glacier length of 3.11 km.The centerline data show that 11 % of the glaciers have advanced in the full epoch (average 88 m), and 5 % and 30 % of glaciers advanced in epochs 1 and 2, respectively.Glacier length changes show strong glacier retreat in all epochs and glacier regions, except for glacier region 19 in the full epoch (Fig. 7b).Overall, the length changes derived in this study show a steady retreat for many of the individual glacier units even though some have advanced (Fig. 7c).It should be emphasized that different numbers of glaciers were compared in the glacier area change and glacier length change assessments because there are significantly less data of the glacier lengths than areas (see Sect. 3.8).This may explain the different patterns of change for epoch 2 (Figs. 5 and 7).Overall, length changes for the three parts of Norway show a retreat of the glaciers for all three epochs (Table 6).

Glacier length changes vs. in situ length changes
Comparison of our glacier length changes with cumulative field data from NVE for 12 selected glaciers shows a mean deviation of 89 m (< 4 %) for the full epoch (Table 4), which indicates that the satellite and map derived glacier length changes are in agreement for groups of glaciers, although for individual glaciers, the deviation can be relatively large.Eight of the glaciers show good agreement (of ±1 to 2 pixels) between the length change methods.
The field measurements reveal that the rate of recession was variable and even absent for some glaciers, varying between −950 m (Fåbergstølsbreen in southern Norway, glacier ID 2289) to +149 m (Engabreen in northern Norway, Glacier ID 1094).Two of the glaciers compared showed a deviation of more than 4 pixels (> 120 m).These glaciers are Fåbergstølsbreen and Stigaholtbreen.
This discrepancy between methods is most likely caused by error in some of the years of the in situ observations (Personal communication with Hallgeir Elvehøy (NVE), December 2013).For all the methods, the local topography in front of a glacier snout can affect glacier length measurements.
Additionally, changes in the glacier's terminus impact the morphology in front of the glacier and make it difficult to compare the measurements (Winkler et al., 2009).The determination of glacier terminus in cast shadow is limited by the quality and resolution of the used satellite images, causing uncertainties in the derived length change (Paul et al., 2011).In our case, Fåbergstølsbreen (Glacier ID 2289) was actually located in cast shadow at the time of acquisition.The deviation can also be caused by differences in the measurement angle of the interpreter on the ground from a reference point toward the glacier and the path of the centerline used to derive length change from glacier outlines.

Glacier change since the beginning of 1900s
Using analogue maps from the beginning of the 1900s (GI 1900 ), we extended the glacier area and length change assessments further back in time for five ice caps or former ice caps (Fig. 8).Similar methods and data sets have been used for detecting a century of glacier retreat on Kilimanjaro (Cullen et al., 2013).This glacier area change analysis included a total of 57 glacier units and 26 centerlines for glacier length analysis, all from these five ice caps.(Table 3).One of the challenges deriving glacier change in this region was disintegrating glaciers, in particular Nordmannsjøkelen (Fig. 8).
The glacier geometry changed extensively, including emerging rock outcrops and ice patches separated from their tributaries.Strong thinning and retreat has been revealed for Langfjordjøkelen, one of the five ice caps, over the period 1966-2008 (Andreassen et al., 2012a).
In total, the five ice caps have decreased in area from 139 to 65 km 2 from the 1900s to 2006 (Fig. 8).The mean area change from GI 1900 to GI 2000 (whole period) for the five ice caps is 1.3 km 2 (53 % in total), and the glacier length measurements show a mean retreat of 1063 m (37 % in total) (Table 3).The mean decadal changes for both area and length show highest retreat rates for epoch 0 between GI 1900 and GI n50 with −0.13 km 2 (length change of −73 m).Epoch 1 and 2 show less relative retreat rates of −0.04 km 2 (length change of −35 m) and −0.03 km 2 (length change of −42 m), respectively (Table 7).Thus, these glaciers have been disintegrating and down-wasting extensively since 1900.

Climate anomalies during the 20th century
Our results show that glaciers in Norway have been receding between GI n50 and GI 2000 , consistent with in situ data of individual glaciers (Andreassen et al., 2005(Andreassen et al., , 2011)).The data also shows that maritime glaciers have been oscillating between periods of advance and recession, with recession being the most frequent state (Andreassen et al., 2005).The strong reduction in area and length in epoch 0 (GI 1900 to GI n50 , Table 7) includes the warm period starting in the 1930s (Hanssen-Bauer, 2005), causing strong glacier retreat (Østrem and Haakensen, 1993;Andreassen et al., 2005Andreassen et al., , 2008)).
Transient mass surplus resulting from increased snow accumulation in the 1990s caused an advance, which is reflected in epoch 2 for both northern and central Norway (Fig. 5a, and see Sect.4.1) as an increase in glacier area (Table 5).The 1990 advance was most likely preserved in the GI 2000 data set in terms of snow around the glacier perimeters and mass gain, especially in maritime areas.
Many of the maritime glaciers have high mass turnover and thus are more sensitive to changes in precipitation rather than temperature (Oerlemans and Reichert, 2000).Glacier variations do not only respond to temperature changes during the ablation season, but are also related to the amount of precipitation in the accumulation season (Nesje, 2005).For several outlet glaciers from Jostedalsbreen, a terminus response time of 3-4 years was observed between the mass gain and the related length change during the 1990s (Winkler et al., 2009).

Elevation
We did not find a significant relationship between average slope over each glacier unit and glacier length or area change for our data sets, although the data show a trend for less glacier change with increasing slope, as previously shown by Leclercq et al. (2014).Our results show that ice caps in northern Norway are particularly vulnerable to glacier area and length changes.Maritime glaciers are in general sensitive in Norway and retreat, but the glaciers in northern Norway retreat more because of less precipitation, warmer temperatures and for many glaciers a location at lower elevations.
(Figs. 6 and 8).These considerable changes are partly attributable to the glacier geometries: ice caps in Norway are relatively flat, and a high fraction of their surface remains close to the modern equilibrium line, which makes them highly sensitive to climatic change (e.g., Nesje et al., 2008), whereas the steep glaciers are less sensitive.If the equilibrium line rises on ice caps, large parts of the accumulation area is transferred to the ablation area, and the mass balance becomes strongly negative.
For example the accumulation-area ratio (AAR) for Langfjordjøkelen, an ice cap in northernmost region, was 0 % for many years during the 2000s, and the glacier is far from being adapted to the present climate conditions (Andreassen et al., 2012a).

Climatic transects
We have inspected how glacier area and length changes are distributed in terms of a west-east transect in our data sets.The west-east transect reaches from Ålfotbreen in the west to Gråsubreen in the east (Fig. 9a).Fig. 9b illustrates yearly glacier area changes in a west-east transect for the full epoch.Annual glacier length changes for the full epoch show similar patterns (Fig. 9c).Even though it is high variability in the data (represented by the boxes in Fig. 9b and c), it shows that the area changes and length changes are consistently more negative for maritime glaciers compared to the continental glaciers.The glaciers located in the precipitation "shadow" of the mountains in the east have less variation in glacier area and length change (Fig. 9).The mean winter balance on Gråsubreen is about 0.8 m water equivalent (w.e.), only ∼ 20 % of the winter balance on Ålfotbreen (3.7 m w.e.) (Andreassen et al., 2005).Representing glacier area changes along a climatic transect illustrates the regional variation of glacier response to climate (Paul et al., 2007).Maritime glaciers have a high mass balance gradient compared to continental glaciers.Due to global warming, glacier retreat will continue, and glaciers in maritime climates are expected to be more sensitive and respond more quickly than glaciers located in continental areas (Hoelzle et al., 2003).The drier continental glaciers are dependent on summer temperatures, while maritime glaciers are more sensitive to spring/fall temperatures (Oerlemans and Reichert, 2000).
Our analysis shows that glacier area and length changes are most pronounced for the northernmost glaciers (Figs. 5 and 7 and Tables 5 and 6).This agrees with geodetic and direct mass balance observations over the last decades.For example, the ice cap Langfjordjøkelen, shows a stronger thinning and retreat than any other observed glacier in mainland Norway.Often the glacier has no accumulation area left at the end of the mass balance year (Andreassen et al., 2012a).The ice cap simply does not have enough area at high altitude for the present climate.
Much of the annual variation in Norwegian climate is influenced by the North Atlantic Oscillation (NAO) (Hurrell, 1995).Glaciers in Norway span a transect of ∼ 1500 km from south to north.Previous studies have shown that the NAO influences the winter and annual surface mass balance, but its effect is reduced towards more continental glaciers, as well as glaciers located at high latitudes (Nesje et al., 2000).

Alternative ways to represent glacier area change
In this paper, we present our data of glacier area changes as absolute (km 2 ), relative (%) and normalized values for different time periods in order to provide a thorough presentation of the data.These methods give different representations of changes.For example, if a small glacier and a large glacier lose the same amount of mass with equal energy inputs, the percent change will be very different between the two.The small glacier will be overestimated in terms of change signal when compared with the large glacier, and the change signal of the small glacier is amplified (Fig. 10b).When glacier change is expressed in terms of km 2 -change, the opposite occurs.A large glacier losing the same area in terms of km 2 as a small glacier with the same energy input, will be overestimated in terms of the change signal, and the signal of small glaciers can be lost (Fig. 10a).Our solution for this was to express glacier area change in a length scale dimension, called units of length.This can be done by normalizing by the square root of the initial area, GAC root(A) , where GAC is glacier area change and A is the initial glacier area which is the GI n50 when representing change for the full epoch (Raup et al., 2009).With this normalization, we removed the systematic trend that depends on initial glacier size (Fig. 10c).By comparing this change signal with glacier length change measurements derived from glacier centerlines, one obtains a similar distribution of the two change signals (Fig. 10d).
Other alternative normalization strategies that also aim at expressing glacier area change in terms of units of length include glacier area change over perimeter or glacier area change over width (personal communication with C. Nuth, May 2013).However, it is important to note that, for Norwegian glaciers, many different glacier types are present, with a variety of sizes and geometries.The "box method" uses the lower part of the glacier tongue for extracting glacier length change (e.g.Moon and Joughin, 2008;Nuth et al., 2013).This method was developed for marine-terminating glaciers, and is not ideal for many Norwegian glaciers, because of the many ice caps and cirque glaciers with often less distinct glacier termini.

Conclusions
We present a glacier area and length change analysis including multi-temporal data sets covering a larger area and higher temporal resolution than earlier studies.The glacierized areas in Norway  a Values are larger than or equal to the given treshold.
b The Blåmannsisen subregion used the threshold values are mapped from three glacier inventories within the period from 1947 to 2006.Overall, glacier area in mainland Norway decreased by −326 km 2 from GI n50 to GI 2000 , corresponding to −11 %.The average glacier length change between GI n50 and GI 2000 is −240 m.Glacier area and length changes indicate that glaciers in western Norway have retreated more than in eastern parts, and northern glaciers have retreated more than southern glaciers.A combination of several factors like glacier geometry and elevation, and climatic aspects such as continentality, are related to the observed spatial trends in the glacier change analysis.The change assessment based on historical maps going back to the 1900s in West-Finnmark revealed a reduction in glacier area from a total of 139 to 65 km 2 , corresponding to a mean area and length change from GI 1900 to GI 2000 of −53 % and −37 %, respectively.Glacier outlines derived from topographic and historical maps have considerable uncertainties due to challenges related to the seasonal snow cover.Therefore, the results show the upper bound of glacier changes in Norway.The results differ regionally, but clearly exhibit a main trend of retreating glaciers between GI n50 and GI 2000 , even though some individual glaciers have advanced.The increased availability of automatically derived and reproducible centerlines makes it easier to retrieve glacier length changes when multi-temporal glacier inventory data are available.Glacier length change derived from centerlines might be a more correct way to express glacier change signals due to reduced dependency on glacier geometries.Sensors with higher spatial and temporal resolution (e.g., the Sentinel-2 satellite) open new possibilities for observing glaciers in the future.ciated in the discussion on alternative ways to represent glacier area change.Many thanks to C. Nuth (UIO)for input and feedbacks on the change analysis in general.Thanks to A. Voksø (NVE) for the design of the geodatabases and rules for topology editing of glacier outlines.Thanks to S. Engh for the manual digitization of the glacier outlines from the topographical maps in his summer job at NVE 2012.Thanks to all colleagues at UIO and NVE for helpful discussions.

Table 3 .Figure 1 .
Figure 1.(a) The study area with Norwegian glaciers shown in blue.(b) Norway is bordered by the Norwegian sea in the west.

Figure 2 .
Figure 2. Spatial representation of the data sets, (a) A subset of five ice caps in northern Norway outlined in the period 1895-1907 (GI1900).The location of the subset is indicated by the black rectangle in 2b.(b) GIn50 consists of 168 N50-map sheets based on aerial photographs within 1947-1985.(c) GI1990 consists of 9 Landsat TM4 and TM5 satellite scenes within 1988-1997.Glaciers not covered by suitable scenes shown in red.(d) GI2000 includes 12 Landsat TM 5 and ETM+ 7 satellite scenes from 1999-2006.(e) Illustration of the division of northern, central, and southern Norway and the 36 glacier regions.

Figure 3 .
Figure 3. (a) Boxplot showing the relative area difference (%) for selected threshold values.The numbers in parentheses are outliers.(b) Visual representation of the threshold sensitivity for a part of Nordmannsjøkelen.The threshold sensitivity is represented in increments of 0.2 using TM 3 TM 5 ≥ 2.4 and TM1 ≥ as the initial threshold values (black line).The blue frame indicates a glacier located in cast shadow.(c) TM TM ≥ 2.4 and TM1 ≥ 60 are used for the initial threshold value.Some of the smallest glaciers or possible snow fields completely vanish when the threshold values are increased.

Figure 4 .
Figure 4. GIn50 and GI2000 presented with the automatic derived glacier heads, termini and centerlines at Øksfjordjøkelen.

Figure 5 .
Figure 5. Glacier area change (GAC) ranging from north to south (a)displayed for Norway, (b)36 glacier regions and (c)1-3143 glacier units.GAC is presented in GAC root(A) (see Sect. 4.5).(a) Boxplot showing the annual GAC root(A) for three parts in Norway, and for three epochs.Only glaciers > 0.5 km 2 are included in (a).(b) mean annual GAC root(A) for 36 glacier regions for the full epoch, and (c) GAC root(A) for each glacier unit for the full epoch.Glacier regions and glacier units are arranged in a north-south order as defined in Andreassen et al. (2012b).

Figure 6 .
Figure 6.Distribution of glacier area with elevation for GIn50 and GI2000 for all glaciers in Norway.We only compare glaciers present in both data sets.The clear bi-modal distribution with a distinction between 1000-1350 m and 1400-1700 m illustrates the predominant location of glaciers in northern and southern Norway, respectively (Andreassen et al., 2012b).Note that northern Norway includes the central and northern part presented in Figs. 5 and 7.

Figure 7 .
Figure 7. Glacier length change (GLC) ranging from north to south (a)displayed for Norway, (b) 36 glacier regions, and (c) 1-3143 glacier units.(a) Boxplot showing the annual GLC (m) for three parts of Norway, and for three epochs, (b) mean annual GLC (m) for 36 glacier regions for the full epoch.Note that glacier regions 16, 18 and 30 do not include any glacier length data, and (c) GLC (m) for each glacier unit for the full epoch.Glacier regions and the glacier units are arranged in a north-south order as defined inAndreassen et al. (2012b).

Figure 8 .
Figure 8. Glacier area and length change for five ice caps in West-Finnmark.The numbers after the mean percentage area and length change are the number of included glacier units (area) or centerlines (length).N = Nordmannsjøkelen, Se = Seilandsjøkelen, L = Langfjordjøkelen, Sv = Svartfjelljøkelen, and Ø = Øksfjordjøkelen.

Figure 9 .
Figure 9. (a) Map of the west-east transect of glaciers in southern Norway, Å = Ålfotbreen, N = Nigardsbreen, and G = Gråsubreen.(b) Mean glacier change every 20 km from the coast for the full epoch is presented in boxplots for yearly glacier area change (units of length GAC root(A) ).(c) glacier length change (m).All glaciers are included for GAC, and glaciers > 1 km 2 are included for GLC.

Figure 10 .
Figure 10.Glacier area change (GAC) and glacier length change (GLC) represented for the Jostedalsbreen region in western Norway.(a) GAC (km 2 a −1 ), (b) GAC (% a −1 ), (c) GAC root(A) a −1 , and (d) GLC (m a −1 ).The corresponding box plots includes all glaciers in Norway > 0.5 km 2 .Extreme values of (a) and (b) are enhanced in the map legend, to mark the point of problematic influence of glacier geometry and size when GAC is represented in km 2 and %.

Table 1 .
The maximum, minimum, and mean time span in years within each epoch.Note the calculated glacier change is weighted by the time span between two data sets for each single glacier.The mean time span in this table is not weighted, but gives the mean of the time span for all glaciers included in each epoch.

Table 2 .
Landsat satellite images for the GI1990 and GI2000 inventories.Dovre, Jotunheimen and Hardangerjøkulen sub-regions were the first sites processed in GI2000, and did not include the TM1 threshold (e.g. Andreassen et al., 2008).(L1G = Image is from the Global Land Cover Facility (GLCF), L1T = Standard terrain correction, NVE = Norwegian Water Resources and Energy Directorate, CGEO = Center for GIS and Earth Observation, ESA (Kiruna) = European Space Agency (Kiruna ground station), USGS = US Geological Survey).