The accelerated decline of mountain glaciers around the world is the clearest natural signal of ongoing atmospheric warming (e.g. Reference Haeberli, Cihlar and BarryHaeberli and others, 2000). According to several climate models, a large temperature increase is expected in particular at high northern latitudes (Reference Moritz, Bitz and SteigMoritz and others, 2002), and corresponding glacier changes are likely to occur. In order to quantify the changes, there is an urgent need for baseline glacier inventory data in these regions, including three-dimensional (3-D) information (e.g. area–elevation distribution for calculation of sea-level rise). Especially high-resolution (10–30 m) multispectral satellite data, such as from the Landsat Thematic Mapper (TM), enable automated and accurate mapping of clean glacier ice over large regions (e.g. Reference Aniya, Sato, Naruse, Skvarca and CasassaAniya and others, 1996; Reference Jacobs, Simms and SimmsJacobs and others, 1997; Reference PaulPaul, 2002a). The possibility of obtaining digital elevation models (DEMs) from sensors with stereo imaging capabilities (e.g. Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER); Système Probatoire pour l’Observation de la Terre (SPOT)) is a most valuable source of topographic information, as the new Shuttle Radar Topography Mission (SRTM) DEM (available at 30–90 m spatial resolution) does not cover these high latitudes (Reference Rabus, Eineder, Roth and BamlerRabus and others, 2003) and contains large data gaps in high-mountain regions. The DEM data are a vital component of any glacier inventory, as they provide 3-D glacier information (e.g. slope, aspect, hypsography).
Within the framework of the Global Land Ice Measurements from Space (GLIMS) project which aims at a global glacier inventory from satellite data (Reference KiefferKieffer and others, 2000; Reference BishopBishop and others, 2004), numerous Landsat 7 and ASTER scenes were provided for analysis to the GLIMS community. In a previous pilot study, the new Swiss Glacier Inventory 2000 (SGI 2000) has revealed algorithms and work flows for efficient glacier inventorying with TM data for mid-latitude mountain glaciers (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002; Reference Paul, Huggel, Kääb and KellenbergerPaul, 2003). Here we apply and evaluate these methods for Arctic glaciers on Cumberland Peninsula, the southernmost part of Baffin Island, Canada (cf. Fig. 1). The focus of this paper is to give an overview of the possibilities for glacier inventorying under Arctic conditions by a synergetic use of satellite and DEM data in combination with Geographic Information System (GIS)-based processing. For a detailed discussion of the respective methods, we refer to the relevant literature. We have also tracked the time required for each processing step, in order to assess the working hours required for a global glacier inventory from satellite data.
In this contribution, we first give a short overview of our study site and then describe the methods and investigations performed. Our focus here is on a comparison of glacier mapping methods, the GIS-based glacier basin delineation, the mapping of Little Ice Age (LIA) maximum glacier extent and the additional capabilities with ASTER data. Afterwards, we present and discuss the results obtained and summarize our conclusions.
Our study site is located on Cumberland Peninsula, Baffin Island, (centred at 66.5˚ N, 65˚ W) with the investigated perimeter defined by a Landsat 7 Enhanced Thematic Mapper Plus (ETM) scene from 13 August 2000 (path 17, row 13) and one adjacent ASTER scene (with small overlap) from the same date (cf. Fig. 1). The region encloses major parts of Penny Ice Cap which is believed to be a remnant of the Laurentide ice sheet (e.g. Reference BirdBird, 1967). For this and other reasons, the area has been the target of numerous field studies during the past 50 years, focusing on glacial– meteorological conditions (Reference OrvigOrvig, 1954; Reference Ward and BairdWard and Baird, 1954), Quaternary history and ice-age chronology (Reference MillerMiller, 1973; Reference DykeDyke, 1979), as well as on equilibrium-line altitude (ELA) and glaciation level (Reference Andrews and MillerAndrews and Miller, 1972; Reference WilliamsWilliams, 1978). The analysis of lake sediments and lichen extent has shed some light on the chronosequence of terminal moraines (e.g. Reference Boulton, Dickson, Nichols, Nichols and ShortBoulton and others, 1976). Moreover, several mass-balance series from Arctic glaciers are available (Reference DowdeswellDowdeswell and others, 1997), and ice-age ELA depressions have been calculated from numerical models (Reference WilliamsWilliams, 1979).
Most parts of the glacierized region are characterized by steep terrain composed of igneous rock with elevations rising from sea level to ∼2300m (Reference BirdBird, 1967). The full spectrum of glacier types is present at our study site (plateau ice field, ice caps, calving glaciers, valley and mountain glaciers, cirques), in part with highly complex shapes and topologies. Although the appearance of numerous glaciers suggests that the ice is inactive (e.g. no crevasses at steep slopes), several glaciers show multiple staggered terminal moraines, indicating highly dynamic glaciers. The location of the terminal moraines with respect to the year 2000 glacier extent is also proof that strong glacier recession has taken place during recent decades in this region (e.g. Reference DowdeswellDowdeswell, 1995). The several push moraines visible are clear evidence for permafrost conditions in the glacier forefield and partly cold glaciers (Reference BennettBennett, 2001). During the 1960s and 1970s, a Canadian glacier inventory was compiled from aerial photography, including Baffin Island (Reference OmmanneyOmmanney, 1980). Detailed analysis of a subsection of Baffin Island can be found in Reference Andrews, Barry and DrapierAndrews and others (1970). Unfortunately, these data are not part of the World Glacier Inventory yet and are not publicly available (Reference OmmanneyOmmanney, 1998). For the above reasons we have selected our study site in this region.
Glacier mapping from multispectral satellite data
Two decades of glacier mapping from multispectral satellite data has resulted in a wealth of methods and data. Although only cloud-free imagery from the end of the ablation period can be used for efficient glacier mapping, appropriate scenes from nearly all the world’s glaciers are available today, due to the long time period covered (starting with TM in 1982). Given that many glacierized regions are not covered by appropriate ASTER/ETM scenes, the archived dataset should also be considered and analyzed systematically for the GLIMS project.
Apart from manual glacier delineation by on-screen cursor tracking, most methods utilize the low reflectance of ice and snow in the middle infrared part of the spectrum for glacier classification (e.g. Reference Bayr, Hall and KovalickBayr and others, 1994; Reference PaulPaul, 2002b; Reference KääbKääb and others, 2003a). This spectral region is covered by several space-borne sensors (e.g. ASTER, IRS-1C/ D, Landsat TM/ETM, SPOT 4/5), and the methods are thus widely applicable. Latest reviews and comparisons of the applied methods were given by Reference Sidjak and WheateSidjak and Wheate (1999), Reference Gao and LiuGao and Liu (2001), Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002) and Reference BishopBishop and others (2004). In particular, thresholded ratio images from various TM band combinations have proven to be a very accurate, fast and robust method, at least for clean glacier ice (Reference AlbertAlbert, 2002; Reference Paul, Huggel, Kääb and KellenbergerPaul and others, 2003). Debris-covered glacier ice is not detectable by multispectral classification alone, but the combination with DEM information (slope) and neighbourhood analysis has shown promising results (Reference Paul, Huggel and KääbPaul and others, 2004 and references therein).
From a remote-sensing perspective, the special challenges at our test site are deep shadows, lakes of varying turbidity adjacent to glacier ice, and debris-covered glaciers. In this section we focus on mapping of glaciers in cast shadow using the ETM scene. In view of the accurate glacier maps obtained from simple band ratios and the absence of an appropriate DEM, we have compared three ratios for glacier mapping: ETM3/ETM5, ETM4/ETM5 and (ETM2-ETM5)/(ETM2 + ETM5), where the latter is also known as the normalized-difference snow index (NDSI). We have also investigated the effect of subtracting path radiance previous to the classification and applied a median filter (noise reduction) to all glacier maps shown in the following comparisons.
In Figure 2a our test region is shown in a contrast-stretched ETM band 3 image (cf. Fig. 1 for location). The most critical regions (shadow, debris, snowfields) are indicated by arrows. The comparisons presented in Figure 2b-d depict the glacier maps from two methods at a time: Figure 2b compares ETM3/ETM5 with ETM4/ETM5, Figure 2c the NDSI with and without subtraction of path radiance in ETM3, and Figure 2d compares ETM3/ETM5 with the NDSI with subtracted path radiance in ETM3. The comparison shows that ETM4/ETM5 is not capable of mapping snow or ice in cast shadow, while the NDSI without path radiance subtraction maps all regions in cast shadow (including bare rock). The visual inspection reveals two ‘winners’ which are compared in Figure 2d, revealing only small differences. The NDSI performs somewhat better in regions with ice in cast shadow, while ETM3/ETM5 is more accurate in regions with snow in cast shadow. Lastly, we use ETM3/ETM5 for glacier mapping, as final manual editing (e.g. for debris cover) has to be applied for many glaciers anyway. Glacier mapping in cast shadow remains the most critical point. Thus, we suggest that efficient glacier mapping benefits most from a careful selection of the threshold parameter in such regions. This will markedly reduce time-consuming manual corrections in these regions.
GIS-based glacier basin delineation
A GIS acts at the interface between the satellite-derived glacier map and the data obtained for individual glaciers. In particular, the digitizing of glacier basins (ice divides) and the automatic data extraction for each glacier, including 3-D parameters obtained from a DEM, are greatly facilitated (Reference Kääb, Paul, Maisch, Hoelzle and HaeberliKääb and others, 2002; Reference PaulPaul, 2002b; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). Here we also use the digitizing capabilities of a GIS to delineate debris-covered glaciers and LIA trimlines or lateral and terminal moraines.
While for the SGI 2000 the digitized Swiss glacier inventory from 1973 was used to define the ice divides (Reference PaulPaul, 2004), such information was not available for Baffin Island. Moreover, the definition of ice divides is hampered, as many small ice caps or larger glaciers are connected in their accumulation area without visible separation. When separation into individual glacier parts was feasible, we used topographic shading to guide the delineation of ice divides. In other regions we used preliminary ice divides (Fig. 3a) as a first step for defining a glacier basin. For a later glacier inventory, these divides should be adjusted according to the former Canadian inventory, which is straightforward with the digital data in a GIS. The divides were extended to closed polygons, roughly surrounding the actual glacier and its possible LIA extent. Contrast-stretched RGB (red, green, blue) composites (bands 3, 2, 1) and the classified glacier outlines are used as background to assist in the digitizing process. With the ETM3/ETM5 method, nearly all lakes are classified as glacier, independent of their turbidity. They are detached from the glaciers in the course of the basin delineation, as they are clearly visible in the background image. Often a glacier lake fills the depression enclosed by the LIA moraine. In this case, the glacier basin polygon encloses the lake, and an additional line is digitized for later lake separation (cf. Fig. 3b). Obvious snowfields within a glacier basin are also enclosed from a polygon and assigned with a label point for later exclusion from glacier area calculation.
Little Ice Age moraines
As a result of the harsh climatic conditions and the recent deposition, the lateral and terminal moraines of many glaciers appear very fresh and are free of any vegetation or even lichen cover. Moreover, most ice caps are surrounded by a lichen-free seam (trimline) of variable extent which is quite obvious even in coarser-resolution satellite imagery (Reference Andrews, Davis and WrightAndrews and others, 1976; Reference Locke and LockeLocke and Locke, 1977). However, when these lichen-free areas are not connected to a recent glacier, they may be attributed to persistent snowfields (Reference KoernerKoerner, 1980). A general recession from the LIA maximum extent did not start until the 1920s, as reported by Reference BirdBird (1967) and Reference DowdeswellDowdeswell (1995). Investigations from lichen dating have shown that the LIA maximum glacier extent (about 130–200 years ago) has been quite similar to mid-Holocene or in some places even late Wisconsinan glacier extent several thousand years ago, and lateral moraines from both periods can be found side by side (Reference MillerMiller, 1973; Reference DykeDyke, 1979). With respect to the sensor resolution of ETM and ASTER (15 m), we are not able to discern these two moraines from space and there might be some interferences in interpretation.
The well-preserved indicators of a former glacier extent form a promising base for manual delineation, in particular as such data have not been obtained thus far for a larger sample of glaciers in this region. However, there are at least three critical points for the delineation:
1. the boundary of the trimlines is quite diffuse (ice caps),
2. there are two or more possible locations for the boundary (lateral/terminal moraines), or
3. they are not visible due to shadow.
In the first case, doubtful parts of the former extent were neglected (Fig. 4a); in the second case, the outline is placed at the inner rim (Fig. 4b); and in the last case, the outline stopped or followed the overall glacier shape (Fig. 4c). Thus, the glacier change calculated here is a minimum value. To delineate former glacier extent, we used pansharpened ETM band 4, 3, 2 RGB composites, obtained by the intensity–hue–saturation (IHS) fusion technique. Here the average intensity from the three 30 m bands 2–4 is replaced with the 15m pan band and than transformed back to the RGB colour space. This composite shows much finer spatial details than without fusion, and the near-infrared (NIR) composite has much better contrast in regions with vegetation than the true colour band 3, 2, 1 composite. The digitizing process includes three steps: the initial digitizing, removal of former misclassified boundaries and correction of dislocated connecting points (node errors) or dangling segments.
Some of the major improvements of the ASTER sensor compared to TM or ETM are: (1) the along-track back-looking NIR band (3B) for generation of DEMs, (2) the improved sensor resolution in the visible/near-infrared (VNIR) bands (15 m), (3) the pointing capability of the VNIR sensor (±24°), and (4) the possibility of changing the gain settings of individual bands (Reference KääbKääb and others, 2003a). Among the disadvantages are the smaller ground area covered (factor of nine compared to Landsat) and the missing band in the blue part of the spectrum (useful for lake detection). Apart from DEM creation, the enhanced capabilities of the ASTER sensor have been utilized in numerous glaciological studies (for an overview see Reference KääbKääb and others, 2003a; Reference BishopBishop and others, 2004): measurement of glacier velocities from sequential images (Reference Dowdeswell and BenhamDowdeswell and Benham, 2003; Reference Kääb, Lefauconnier and MelvoldKääb and others, 2005), assessment of glacier hazards (Reference Huggel, Kääb, Haeberli, Teysseire and PaulHuggel and others, 2002; Reference Kääb, Wessels, Haeberli, Huggel, Kargel and KhalsaKääb and others, 2003b), analysis of supraglacial lakes (Reference Wessels, Kargel and KiefferWessels and others, 2002) and mapping of debris-covered glaciers (Reference Paul, Huggel and KääbPaul and others, 2004). The generation and accuracy assessment of ASTER-derived DEMs has been discussed in detail before (e.g. Reference KääbKääb, 2002; Reference ToutinToutin, 2002a, b) and is thus only briefly summarized below.
We apply ASTER data in this study for two major objectives: generation of a DEM and comparison of the glacier mapping results with ETM. While glacier mapping within the ASTER scene is somewhat hampered by snowfields, the DEM creation benefits from the low gain setting in the VNIR bands. The DEM has several potential applications, including orthorectification of satellite imagery, 3-D glacier parameters, perspective visualizations, topographic normalization, potential solar radiation, identification of ice divides, and mapping of debris-covered glaciers. Here we apply the DEM for orthorectification and calculation of 3-D glacier parameters for 340 glaciers with the method described in Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002). The semi-automatic method for delineating debris-covered glaciers (cf. Reference Paul, Huggel and KääbPaul and others, 2004) did not work well in this region and has thus not been applied. For the ETM scene (no DEM available) we applied manual delineation of debris cover for 85 glaciers using a pan-sharpened band 4, 3, 2 composite as background.
The DEM was obtained from the ASTER 3N and 3B channels at 30 m spatial resolution using 17 ground-control points (GCPs) and 25 tie points within standard image-processing software (PCIs Orthoengine). Elevations and locations of the GCPs are taken from the 1 : 250 000 scale Canadian topographic map (sheet 16L) for several mountain peaks and some points at sea level. The average difference of calculated and mapped elevations is –8.6m, with a standard deviation of 155 m. Visual inspection of the shaded DEM reveals only few locations with outliers, mostly located at steep, shadowed slopes and not on glaciers. Small gaps in the DEM are filled by interpolation in the course of geocoding. High-frequency noise in the resulting DEM has been reduced by a Gaussian low-pass filter which changed the original elevations by only a few metres.
Glacier mapping with ASTER
The glacier map from ASTER has also been obtained from a simple band ratio using a thresholded band 3/4 image and an additional threshold in band 1. The latter removes nearly all of the wrongly classified shadow pixels that are located over bare rock, due to the markedly higher reflection of shadowed ice and snow pixels. For a small section shared with the ETM scene, we compared the glacier mapping results of both sensors. Visual inspection reveals very good agreement of the two outlines (Fig. 5a), apart from small shifts due to the missing orthorectification of the ETM scene. A more quantitative evaluation has been performed for 62 glaciers >0.1 km2, yielding an average relative area difference of -2.4%, with a standard deviation of 5.6% (Fig. 5b). Thus, ASTER-derived glacier areas are in general somewhat larger, but the deviations can also be explained by differing pixel sizes (Reference Paul, Huggel, Kääb and KellenbergerPaul and others, 2003). For glaciers larger than 1 km2 the differences are within ±5%. For clean glacier ice, previous studies have revealed area mapping errors with TM of <5% (e.g. Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others 2002, Reference Paul, Huggel, Kääb and Kellenberger2003). However, this error is largely theoretical, as nearly every glacier is partly covered by debris.
The distribution of glaciers by number and area for seven distinct area classes is depicted in Figure 6a for 770 glaciers as obtained from the ETM scene. While the percentage of glacier-covered area clearly increases towards larger glacier size classes, the percentage of glaciers by number is more-or-less similar for all size classes less and larger than 5 km2, respectively. For Alpine glaciers, nearly the opposite is the case, with more-or-less similar percentages by area for all size classes (summing up all glaciers smaller than 1 km2) and a nearly exponential decrease of the percentage by number towards larger glaciers (Reference PaulPaul, 2004). Thus, we suggest that calculation of global glacier distribution from universal scaling paradigms (Reference BahrBahr, 1997) may give erroneous results. There is not much change in this general pattern if size classes are grouped in powers of two as in the World Glacier Inventory.
Area change since the LIA maximum extent
The relative change in glacier size with area is illustrated for 225 individual glaciers in Figure 6b. The relative loss of area increases towards smaller glaciers (if sizes are <10 km2) and is more-or-less constant for larger glaciers. Interestingly, there is no increase in scatter towards smaller glaciers as observed, for example, in the Swiss Alps (Reference PaulPaul, 2004). In total, the 225 glaciers have lost 217 km2 of area (–11%) since their LIA maximum extent. As several authors have pointed out, general glacier recession in this area started in the 1920s (e.g. Reference DowdeswellDowdeswell, 1995), resulting in an average decadal area loss of 27km2, or 1.4%. This value is somewhat lower than the 38 km2 (–2.2%) area loss per decade found for Alpine glaciers for the period 1850–1973, but still indicates that glacier recession in this region is considerable.
Several 3-D parameters were obtained for 340 glaciers from the ASTER-derived DEM and the glacier outlines obtained according to the method described in Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002). In order to keep this overview short, the two scatter plots in Figure 7 are presented only as an example of this dataset. Future studies might explore the full potential of the DEM for further applications. From Figure 7a a clear dependence of minimum and maximum glacier elevation on glacier size is evident, although there is some scatter towards smaller glaciers. The scatter plot in Figure 7b illustrates about 300m higher mean elevations for south-oriented than for north-oriented glaciers. For all 340 glaciers, the average of the mean elevation is 1020 m (±217 m), the elevation range is 550 m (±385 m) and the mean slope is 23° (±7.7°).
In view of the GLIMS goal to create a global glacier inventory from space, we have tracked the time required for glacier mapping, in particular glacier basin delineation and manual editing of glacier outlines (debris cover, shadow gaps). We have spent about 22 hours on basin delineation for 880 glaciers (1.5 min per glacier) and about 4 hours on debris mapping and final editing of 85 glaciers (3 min per glacier). This sums to an average 5 min per glacier, or 7 years for all estimated 160 000 glaciers worldwide (based on 240 8 hour days per year). Mapping of LIA outlines took an additional 20 hours for 230 glaciers (5 min per glacier). The mapping time required for a full scene with the ratio method is about 1 hour, including algorithm and threshold selection for one or two test regions. Thus, glacier mapping is very fast, but post-processing takes considerable time.
While glacier basin delineation is straightforward for most cirques, as well as for valley glaciers which are separated by mountain crests (as in the Alps), there are problems defining the ice divides of many glaciers found in the Arctic. For assessment of glacier change, we have introduced preliminary divides which can be corrected in a later revision. LIA glacier outlines have been delineated because they are prominent in the satellite imagery and it is valuable to quantify glacier retreat in this area. For some LIA glacier outlines, we cannot exclude misinterpretation in this first assessment, but the overall pattern and magnitude of glacier retreat will not change significantly after a later refinement of glacier basins. The 15m ASTER VNIR bands show finer spatial details (cf. Fig. 5a) than the corresponding pansharpened 15m ETM bands. This permits improved delineation of LIA moraines and debris-covered glaciers. The numerous applications of the ASTER-derived DEM were not fully exploited in this overview. However, the examples shown illustrate the potential of such DEMs, even if the cartographic database is coarse (1 :250 000 scale map). From a GLIMS perspective, the ASTER DEM is a vital part of the project, and the combination with the SRTM-derived DEM will be most fruitful in other parts of the world.
Conclusion and Perspectives
We have discussed the methods and possibilities for glacier inventorying from space by the synergetic use of ETM and ASTER data at a remote site in the Canadian Arctic. While accurate glacier mapping can be obtained quite fast from simple band ratios, the amount of work involved in postprocessing (delineation of glacier basins, debris-covered parts and misclassification) should be taken into account (5 min per glacier). The ASTER-derived DEM in combination with GIS-based processing provides a wealth of additional possibilities, including automated 3-D glacier inventorying, particularly for regions outside the coverage of the SRTM DEM. In future studies we will extend the sample presented here to adjacent ASTER as well as ETM scenes and hope to integrate the data obtained in the GLIMS database at the US National Snow and Ice Data Center.
All ASTER and ETM scenes used in this study are kindly provided by the GLIMS science team. We would also like to thank T. Kellenberger for maintaining the latest image-processing software at our institute. We gratefully acknowledge the helpful comments of the anonymous referees and the suggestions made by the scientific editor J. Dowdeswell.