Trimline Mapping from Multispectral Landsat ETM+ Imagery

Csatho, Bea; Van Der Veen, Cornelis; Tremper, Catherine

doi:https://doi.org/10.7202/013736ar

Introduction

Constructing longer glacier histories, spanning the last few hundred years or so, is of particular interest for outlet glaciers draining the interior of the Greenland Ice Sheet. Repeat surveys by airborne laser altimetry in the 1990s have revealed significant thinning of the ice sheet at lower elevations (Krabill et al., 2000). Specifically, some of the larger outlet glaciers are thinning at rates of several meters per year, with the most extreme thinning observed near the terminus of the Kangerlussuaq Glacier in southeastern Greenland and over the lower reaches of the Jakobshavn Glacier (Jakobshavn Isbrae) in west Greenland, with measured thinning rates as high as 10 m/yr in both locations (Thomas et al., 2000, 2003). Such large thinning rates cannot be explained by a recent increase in ablation or decrease in snowfall and, instead, point to ice-dynamical processes as drivers of glacier retreat (Van der Veen, 2001; Mitchell, 2005). To fully appreciate the significance of these recent glacier changes, the magnitude of retreat and surface lowering must be placed within the broader context of retreat since the Last Glacial Maximum (LGM) and, more significantly, retreat following the temporary glacier advance during the Little Ice Age (LIA). The instrumental record of glacier observations in Greenland dates back to aerial photography conducted by the Danish Geodetic Institute in the 1930s and 1940s. Glacier histories extending farther back in time must be based on geological information retrieved from formerly glaciated regions. Depending on the age of mapped glacial-geological features, histories may be extended to the LIA and possibly as far back as the LGM, some 18 000‑21 000 years ago.

The traditional approach to reconstructing paleo ice sheets from geomorphological and geological information has been to map such features during detailed field investigations. Clark (1997) noted several barriers inherent to this approach that impede synthesizing such fragmentary information into a coherent ice-sheet wide picture. In most cases there are few links between spatially separate studies, or there are large gaps in spatial coverage. Moreover, in many instances, control on dating is poor, thus casting doubt on the commonly-made assumption that similar features reported in different studies are coeval. These difficulties can be partially overcome through the use of remote sensing, which allows for greater spatial coverage of geological features (Clark, 1997).

As he has noted, “A major problem in validating ice sheet models is that geologically-based information has been derived mostly from fieldwork and usually relates to only small sectors of ice sheets which can be thought of as point evidence. Furthermore there are frequently unresolved contradictions between local areas. A methodology that attempts to provide information on wider spatial scales and promote greater analysis and coherence of all evidence is by use of remote sensing and … GIS”.

Early geological applications of remote sensing were based on high-resolution aerial photogrammetry. According to Clark (1997), the study of Sugden (1978) was the first to demonstrate the usefulness of satellite images to provide glacial-geomorphological information over large areas. In that study, a Landsat MSS photomosaic was used to produce a proxy map of the intensity of glacial erosion by the Laurentide Ice Sheet. Subsequent applications have included mapping of depositional landforms (Punkari, 1982), identification and mapping of previously unrecognized mega-lineaments (Clark, 1993; Clark et al., 2000), and mapping margins of paleo ice streams (Stokes and Clark, 2002). A review of earlier work is provided by Clark (1997). The objective of the present study is to explore the feasibility of using multispectral Landsat ETM+ images to map trimlines in the coastal regions of Greenland.

In mountainous regions, former higher stands of glaciers that retreated recently (since the Little Ice Age maximum) are frequently marked by trimlines on the walls of valleys and fjords previously occupied by these glaciers (Fig. 1). Erosional action of the moving ice causes vegetation to be removed where the ice is in contact with the surrounding rock walls. The upper limit of this eroded region, which becomes exposed as the glacier retreats and its surface is lowered, is called the trimline (Flint, 1971: p. 144). Mapping of trimlines and other glacio-geomorphological features thus allows reconstruction of the neoglacial history of glaciers.

Historical evidence indicates that in many parts of the Greenland Ice Sheet the major historical maximum occurred almost simultaneously during the last decades of the 19^th century (Weidick, 1984). Therefore, mapping of the trimline zone over larger areas and possibly the entire ice sheet could provide an estimate of the total area deglaciated as well as the total ice volume lost since then. Indeed, based on comparison between trimline elevations and the current ice-surface elevation in west Greenland, Weidick (1968) estimated an average loss of 200 km³/year ice from the whole Greenland Ice Sheet since the Little Ice Age, corresponding to an average rise of 0.6 mm/year in global sea level. This is equivalent to a 30 800 km³ total loss or 90 mm global sea level rise since 1850.

The ability to map surficial features from satellite imagery depends on obtaining spectral signatures of these features that are markedly different from the signatures of surrounding regions. With respect to trimlines, spectral reflectance differences may be associated with variations in vegetative cover above and below the trimline. The LIA trimlines are usually characterized by sharp vegetation boundaries. Outside the trimzone, flat and gentle slopes are often covered by tundra vegetation, while steeper rock surfaces are lichen-covered. Inside trimline regions, rock surfaces are fresh and bare, retaining their original color and spectral properties (Knight et al., 1987). It is well known that vegetation can easily be mapped using visible and near infrared (NIR) spectral imagery. Several studies have demonstrated the pronounced effect of lichens on the reflectance spectrum of the land surface (Satterwhite et al., 1985; Ager and Milton, 1987; Rivard and Arvidson, 1992; Bechtel et al., 2002; Van der Veen and Csatho, this issue). Thus, spectral property information from multispectral satellite images can be used to locate the boundary between vegetation-covered and vegetation-free surfaces. This boundary is the trimline.

Characteristics of lichens

Lichens are “an association of a fungus and a photosynthetic symbiont resulting in a stable thallus of specific structure” (Hawksworth and Hill, 1984). Lichens occur in most habitats around the world and communities are usually classified according to their substrate. Worldwide, approximately 13 500 taxa of lichenized fungi are recognized (Hawksworth and Hill, 1984). Lichens, which do not have special organs or structures for water storage, shut down metabolically when dry, allowing them to withstand longer periods of dryness. Because of this opportunistic usage of available moisture, lichens can survive in environments that are low or deficient in moisture available for the support of plant life. In these so-called xeric environments, lichens spend much of the time in a dry, inactive state (Lawrey, 1984).

Rock substrates provide stable habitats for the establishment of long-term lichen communities. On these surfaces, lichens grow at relatively uniform rates of time scales of centuries or millenia. Once attached to a substrate, lichens do not move during their lifespan and thus the age of lichens can be considered a proxy for the minimum exposure time of the substrate to the atmosphere and sunlight. These two characteristics form the basis for lichenometry, a biological technique for estimating the age of exposed and stabilized rock surfaces (Locke et al., 1979). In short, the assumption is made that the size of the largest individual thallus on a substrate is a function of the time the substrate has been exposed. Using surfaces of known exposure age, a lichenometric dating curve can be established that relates surface exposure age to the largest lichen size. This dating curve is then applied to estimate the exposure age of other surfaces based on the largest size of lichens found on these surfaces (Noller and Locke, 2000). Lichenometry has been applied in glacial geology and geomorphology at many sites around the world including dating the LIA maximum in Iceland (Kirkbridge and Dugmore, 2001), in South America (Winchester and Harrison, 2000) and in Alaska (Haworth et al., 1986). Beschel, who is regarded as the “Father of lichenometry”, has pioneered lichenometric dating in west Greenland (Beschel, 1961).

While the study presented here focuses on determining the spatial extent of lichen-covered rocks, exploratory lichenometry was conducted and a tentative dating curve for dating glacio-geomorphological features exposed since the LIA was derived. Results of that part of our study will be discussed in a forthcoming paper.

The spectral characteristics of lichens are discussed more fully in Bechtel et al. (2002) and in a companion paper (Van der Veen and Csatho, this issue), but some details are worth repeating here. Spectra of lichens in the 400‑700 nm range vary primarily with pigmentation, and in the 700‑1 350 nm range reflectance increases steadily and remains high or only slightly decreasing at greater wavelengths, in contrast to the reflectance of green, leafy plants, which tends to rise almost vertically from the visible to the near-infrared and decreases slowly from 800‑1 300 nm (Ager and Milton, 1987). Absorption bands centered at 1 450, 1 935 and 2600 nm are attributed to water in lichens. Since these water absorption features are usually weaker than in vascular plants, they do not mask the effect of cellulose, which causes broad absorption bands at 1 730, 2 100 and 2 300 nm (Bechtel et al., 2002).

To our knowledge, the only previous study of spectral properties of lichens and lichen-covered rock from Greenland was published by Rivard and Arvidson (1992). The authors collected in situ spectra of different surfaces, including lichen- covered rocks, on Qilangarssuit Island on the west coast to investigate whether Landsat Thematic Mapper (Landsat TM) multispectral imagery can be used for lithologic mapping in Arctic regions. They concluded that Landsat TM data can be used as first order products to distinguish between tundra vegetation that typically covers older moraines at lower elevations, and lichen-covered bedrock exposed at higher elevations. Mapping of different lithologies with Landsat TM imagery is severely hampered by the limited spectral and radiometric capabilities of the sensor as well as by the ubiquitous lichen cover. For the purpose of mapping trimlines, however, requirements are less stringent and the results of Rivard and Arvidson (1992) are encouraging as these indicate the possibility of distinguishing between tundra cover, lichen-covered rocks, and bare rock surfaces.

Previous trimline studies

It has long been recognized that trimlines can be mapped using aerial photography. For example, Weidick (1992; reproduced in Fig. 7: right panel) presents glacial erosion trimzones along the margins of Jakobshavn Isbrae in west Greenland, mapped from stereo aerial photographs acquired in 1985. However, as noted by Clark (1997), aerial photographs, while capable of providing excellent high-resolution information on geomorphological features, typically cover limited areas and extending this procedure to cover regions formerly occupied by entire ice sheets is a daunting task. Because of the greater area captured by a single image (typically 60 by 60 km or larger), satellite imagery permits greater areas to be studied efficiently. To our knowledge, the study of Knight et al. (1987) is thus far the only study to use satellite images to map trimlines and ice-sheet retreat.

Landsat Multispectral Scanner (MSS) imagery, acquired in July 1982, was used by Knight et al. (1987) to identify recently deglaciated regions bordering the present-day ice sheet margin in the vicinity of Jakobshavn Isbrae. In the MSS-5 band (red, 600‑700 nm), a light-tinged zone fringing the ice sheet and extending down both flanks of Jakobshavn Isfjord was noted. The outer limit of this zone appeared to be marked by a sharp boundary in vegetation cover. Using supervised classification of all four MSS bands, Knight et al. (1987) distinguished six different classes of surface types, namely sediment filled lakes, clear lakes, ice cap, glacier debris, vegetated area, and poorly vegetated areas within the trimline. Feature space boundaries of the classes for the parallelepiped classification were defined by the standard deviation of the spectral values in the training sets. The good separability of the classes in feature space indicated that sufficient spectral contrast exists between fresh rock surfaces, and vegetation and lichen-covered rocks. In particular, poorly vegetated areas within the trimline were isolated from regions covered with shrub tundra and lichen cover, primarily through their higher reflectance in bands MSS-4 (green, 500‑600 nm) and MSS-5 (red, 600‑700 nm). Since the spectral contrast along the trimline depends on the type of vegetation, rocks and moraines exposed within the trimzone (Satterwhite et al., 1985; Ager and Milton, 1987; Van der Veen and Csatho, this issue) the four bands of the Landsat MSS sensor allow detection of trimlines only for certain combinations of lichen and rock reflectance. Therefore, the procedure followed by Knight et al. (1987) may not be directly applicable to other sites. The more recent Landsat TM and Enhanced Thematic Mapper (ETM+) systems have six spectral bands within a wider spectral range (450 to 2 350 nm) enabling a better discrimination between different landcover types and therefore these sensors are more suitable for mapping trimlines. Moreover, the Landsat TM and ETM+ systems have higher spatial resolution, namely 30 m Instantaneous Field Of View (IFOV) of the multispectral bands and 15 m IFOV of panchromatic band (ETM+ only) compared with the 79 m IFOV of the Landsat MSS, enabling mapping of smaller geomorphological features, such as small moraines.

Feasibility study: Jakobshavn Isbrae

The drainage basin of Jakobshavn Isbrae (known also as Sermeq Kujalleq), named after the nearby town of Jakobshavn (Ilulissat) on the west coast of Greenland (Fig. 2), is one of the most dynamically active parts of the Greenland Ice Sheet. The calving terminus of this glacier receded 25 km between 1850 and 1950 from the head of the fjord to the position shown on the Landsat image used in this study (Weidick, 1992: Fig. 2; reproduced in Fig. 7). Fresh moraines and trimlines exposed on the fjord walls suggest a lowering of the glacier surface by up to 250‑300 m since the LIA maximum (Weidick, 1992). Recent results indicate further retreat (Podlech and Weidick, 2004), widespread thinning over the lower part of the drainage basin (Thomas et al., 2003) and substantial increase in velocity (Joughin et al., 2004). Because of these ongoing changes, as well as the availability of prior studies in this region (Thomsen, 1983; Thomsen, 1986; Knight et al., 1987; Thomsen and Winding, 1988; Thomsen et al., 1988; Weidick, 1992; Prescott, 1995; Hughes, 1998; Sohn et al., 1998), this area was selected for the feasibility study.

The Landsat ETM+ imagery, taken on July 7, 2001, covers Jakobshavn Isfjord and a 200 km long stretch of the ice sheet margin north of the fjord (Fig. 2). The boxed area around the fjord, measuring approximately 65 km by 75 km, was selected for this pilot study because of the availability of detailed land cover information. In addition to image analysis described in this paper, exploratory fieldwork was conducted around Ilulissat and in and around the field camps marked in Figure 3. The major objectives of the field campaign were to collect the spectra of various land surfaces and vegetation types and to map glacio-geological features, including trimline locations. Van der Veen and Csatho (this issue) summarizes the results of the spectral measurements. Results from the spectral measurements and from the glacial geologic mapping were also used in this study, for example for training site selection and for assessing the accuracy of the classification. The study site is depicted in Figure 3 as a false color composite of Landsat ETM+ bands 2, 3 and 4. This “color infrared” image shows vegetation in different hues of red, while water, snow, and ice surfaces are depicted in blue, ranging from the dark blue over sea water to white over snow covered regions. Exposed rock surfaces and man-made objects are characterized by grayish-brownish hues. At the time of the image acquisition, most of the seasonal snow cover had melted from the coastal mountains, exposing the rock and vegetation surfaces.

The land area is hilly upland, composed of rounded knolls of crystalline bedrock with elevation ranging from sea level to ~650 m. The bedrock consists of Precambrian crystalline rocks, mostly gneiss, with amphibolite dykes in the immediate study region (Geological Survey of Denmark and Greenland, 1997). The area, typical of west Greenland, is characterized by a landscape of glacial scouring (Sugden, 1974). The satellite image depicts the landforms of glacial erosion (Fig. 3), irregular depressions having been produced along joints, faults and dykes, some of which are occupied by lakes in the deeper basins.

Over the exposed land the predominant reflectance feature is the high NIR reflectance (reddish color of the IR composite) indicating that the surface is covered by a patchwork of vegetation and lichens. The vegetation is of low-arctic type with predominant plant communities of heath (moor), fell-field, herb-slope, willow-scrub, fen, river- and seashore vegetation. Mikkelsen and Ingerslev (2002) include details on these plant communities as well as a comprehensive plant list of the Ilulissat area. The most abundant lichens are dark grey to black crustose and foliose ephilitic (rock-growing) lichen (Van der Veen and Csatho, this issue). The brownish-grayish band between the ice sheet and the land corresponds to fresh rock surfaces exposed by retreat of the ice sheet since the LIA. Similar light-colored bands are observed around some of the lakes, where the former higher water level destroyed the lichen cover.

The clear water in fjords bordering Disko Bay (Disko Bugt), and supraglacial lakes are shown in dark blue. Most of the lakes in the ice free terrain also have dark blue colors. Lakes close to the ice sheet margin have lighter blue shades indicating higher reflectance in the visible and in the NIR, probably caused by fine-grained, suspended sediment, originating from englacial and subglacial channels and glacier melt-streams. Fresh snow covering the higher elevation parts of the land and ice sheet surface is white. Older snow, firn and glacier ice appear in different shades of blue. The image also depicts the details of the ice sheet surface morphology, including crevasse zones, meltwater streams, and debris cover adjacent to the margin. Jakobshavn Isfjord is filled with brash ice surrounding a few icebergs and bergy bits. Some of the lakes on the surrounding land are still covered by ice.

Analysis of Landsat imagery

Preprocessing

Orthorectification

The orthorectified Landsat ETM+ scene, taken on July 7, 2001 (Fig. 2), was obtained free of charge from the Global Land Cover Facility (http://www.landcover.org). The sun elevation angle was 42 degrees and the sun azimuth was 169 degrees during the image acquisition. Landsat images of Greenland in this global data set were registered by using manually selected control points in a block-adjustment involving all scenes covering Greenland, followed by orthorectification using the 30‑arc-second GTOPO30 (1 km) Digital Elevation Model (Tucker et al., 2004). The estimated planimetric accuracy of the orthorectified images is better than 50 m Root Mean Square Error (RMSE).

Atmospheric correction

We applied the Empirical Line Correction (EML) model to convert the raw Digital Number (DN) values of the six multispectral bands into surface reflectance. The EML technique uses a band-to-band linear regression between the pixel samples and reference reflectance spectra to establish the transformation from the DNs to surface reflectance (Kruse et al., 1990). Because of the remoteness of sufficiently large homogeneous regions, we could not collect field spectra suitable for establishing this transformation. Therefore we followed the normal practice of using laboratory spectra as reference spectra (Rosenthal and Dozier, 1996) by selecting water, snow and gneiss laboratory spectra from the Johns Hopkins University Spectral Library (included in Environment for Visualizing Images, ENVI 4.1, from Research Systems, Inc.). Corresponding regions of the Landsat imagery are marked in Figure 3. Water pixels were extracted over open sea water. For snow, we selected an area of darker, possibly melting, snow to avoid saturation of the visible bands over bright, fresh snow. These snow pixels, taken near the ice margin, SE of Jakobshavn fjord, were paired with the spectra of coarse, granular snow from the spectral library. Landsat pixels over outcropping rocks, south of Sermeq Avannarleq, and felsitic gneiss reference spectra were chosen for representing the typical granodioritic to quartz dioritic gneiss rocks of the region (Geological Survey of Denmark and Greenland, 1997). To assess the performance of the EML conversion we compared the reflectance values derived from Landsat imagery with spectra measured by using a Fieldspec^® PRO spectrometer from Analytical Spectral Devices, Inc. over the runway of the Ilulissat airport during our 2003 field campaign. The good agreement between the reflectance values [reflectance_Landsat = (0.089, 0.117, 0.123, 0.145, 0.190, 0.194) and reflectance_field = (0.110, 0.125, 0.129, 0.139, 0.175, 0.185), vector components refer to Landsat spectral bands] confirms the accuracy of the EML reflectance conversion.

Topographic correction

A Digital Elevation Model (DEM) of the site, generated from stereo Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) by using the ASTERDTM ENVI module from Sulsoft, Brazil, was used for topographic correction. Stereo aerial photographs and GPS observations, both from the National Survey and Cadastre, Denmark, provided control information for the generation of the DEM. The application of a simple topographic correction procedure based on the cosine law (Bishop and Shroder, 2004) caused overcorrection over steep slopes. This effect is especially pronounced over slopes with cast shadows facing away from the sun. This problem occurs because the cosine correction does not account for secondary illumination, such as diffuse skylight and light reflected from adjacent terrain. Since the use of the topographic correction did not increase the classification accuracy, as was demonstrated by several test runs over small sites, we decided to use the uncorrected reflectance images.

Surface classification

Before the supervised classification of the non-thermal Landsat bands, we examined their spectral information content by using unsupervised classification and different visualization tools (e.g. scatterplots and n-dimensional visualizer) provided by the ENVI software. We identified 14 different spectral classes, related to different snow, ice and water conditions, and different vegetation, rock and sediment types (Table I). Training sites are well distributed within the study area shown in Figure 3. When selecting the training sites we relied on the following data sources: (1) topographic maps (National Survey and Cadastre of Denmark, 2000); (2) geomorphologic maps (Thomsen and Winding, 1988); (3) geologic maps (Geologic Survey of Greenland, 1974; Geologic Survey of Denmark and Greenland, 1997); (4) stereo aerial photographs acquired in 1985 (National Survey and Cadastre of Denmark, 2000; 1 : 150 000 scale); (5) historical photographs, for example in Mikkelson and Ingerslev (2002) and in Weidick (1968); (6) field observations conducted around the town of Ilulissat, around Camps 1‑3 and in the Pâkitsoq area in 2003 (Fig. 3).

Separability of classes was checked during the selection process by computing the transformed divergence of the training sets (Schowengerdt, 1997). Transformed divergence measures the separabilty of samples by considering both their means and covariance matrices. If the transformed divergence of classes was less than the threshold of good separability (1.9), the training sets were modified. A total of ~108 200 training pixels in 59 training areas (~2% of the pixels in the classified scene) were selected to establish the statistics of the 14 landcover classes.

Table I

**Major spectral and information classes of the study site**

The classification map, presented in the upper panel of Figure 4, was obtained by Maximum Likelihood (ML) classification of the six non-thermal multispectral bands of the Landsat ETM+ imagery. The classified image was filtered by using a 3 by 3 majority filter (weight of central pixel = 3), to remove isolated pixels. To delineate the trimline zones, the spectral classes were combined into five information classes in the lower panel of Figure 4: (a) ice, snow, firn and supraglacial lakes; (b) water (from clear to turbid); (c) land, including vegetation-covered surfaces and glaciofluvial and eolian deposits; (d) debris-covered ice and (e) trimzones.

The average spectra of the 14 different surface classes are depicted in Figure 5. Although the low spectral resolution of the Landsat sensor and the effect of spectral mixing somewhat complicate the interpretation of these spectra, most land cover type exhibits well-known spectral shapes and absorption features.

Freshly fallen snow has a very high reflectance in the visible wavelength and the reflectance rapidly decreases toward the NIR and MIR (spectrum A). Spectrum B, observed mostly over upland lakes and in fjords, has a shape similar to spectrum A, but the reflectance is lower throughout the whole spectral range. Reflectance of fresh snow, firn, and ice is reduced by the presence of impurities, debris cover and surface meltwater (Hall and Martinec, 1985), and we interpret class B as melting ice or pieces of floating, calving ice in fjords. The shape of spectrum C, which includes pixels measured over supraglacial lakes, is quite different. The high reflectance at blue and green (bands 1‑2) indicates that the energy at shorter wavelengths penetrates deep enough to be reflected from the bottom ice. However, as the penetration depth in water decreases with increasing wavelength, the reflectance is rapidly decreasing and the spectrum becomes similar to water spectra.

Open, clear water (spectrum D) has a low overall reflectance. Multiple scattering, caused by suspended silt in the water, increases the reflectance in the visible and NIR wavelengths observed in spectra E, F and G, with increasingly higher reflectance indicating progressively more turbidity (Mertes et al., 2004: p. 374). The ability to distinguish changes in turbidity could permit monitoring sediment load in proglacial lakes, streams and fjords from multispectral satellite imagery. Since most suspended sediments originate from englacial and subglacial channels and glacier melt-streams, this, in turn, could provide important information about glacial and subglacial hydrologic systems.

Spectra of biogenic surfaces (spectra H, I and N) show typical spectral signatures of lichens (the reflectance maximum in the MIR band) and vascular vegetation (sudden increase of reflectance in NIR, also called “red edge”). Class I shows a more pronounced red edge, with a smaller maximum in band 5, suggesting that this class has more tundra vegetation and less lichen-covered rock than class H. To explain the spectral signatures of these classes, we computed the spectral curves of mixed pixels from field spectra of typical land cover types, collected around Camp-2 during the 2003 field season (see Van der Veen and Csatho (this issue) for details in data acquisition). The upper left panel of Figure 6 shows the field spectra of (1) Slip arctica (Arctic willow) that represents vascular vegetation; (2) Umbilicaria cylindrica, a black, foliose lichen, covering large parts of rock outcrops above the trimline and (3) light-colored gneiss, which is a significant constituent of the exposed soils and sediments. The spectral library resampling tool of ENVI 4.1 was used to resample these field spectra into the spectral resolution of the Landsat ETM+ sensor, resulting in the spectra presented in the upper right panel. We computed the spectra of mixed pixels, shown in the lower panel (thin lines), by assuming linear mixing of the different components. By comparing these spectra with the average spectra of vegetation-covered areas from the Landsat ML classification (shown by thick lines), we estimate that class H (lichen-dominated vegetation) has 70% lichen and 20% vascular vegetation cover, while class I (tundra and lichen vegetation) has only slightly more lichen-covered area (50%) than vascular vegetation (40%). The percentage of outcrop is 10% for all modeled cases. Our results demonstrate that Landsat TM and ETM+ data can be used as first order products to distinguish between tundra vegetation that typically covers older moraines at lower elevations (class I), and lichen-covered bedrock exposed at higher elevations (class H). We should emphasize that determining the distribution of different land cover types within mixed pixels is an ill-posed problem and therefore the result of spectral unmixing is somewhat influenced by the end-member selection.

The spectra of exposed sediments and rocks (spectra J, K and M) have very similar shapes with no major absorption features. Class K, characterized by a reflectance increasing with wavelength in the visible and near constant, high reflectance in the NIR/MIR domains indicates the presence of fine grained material, with low organic matter content and with no or minimum vegetation cover. This class is found in outwash plains and fluvial terraces deposited from braided rivers and according to the Quaternary geologic map it corresponds to eolian and glaciofluvial sediments. The spectrum of class J has a similar shape but an overall lower reflectance with a pronounced decrease in reflectance toward the end of the spectrum, indicating similar material but higher moisture content. Therefore class J is interpreted as wet glaciofluvial and eolian sediments. The featureless spectral curve of M is typical for the moraines and outcrops of the trimline zones. Finally, the spectrum of class L is similar to the spectra of moraines and the trimline zone (class M) in the visible and NIR domains, but has slightly lower reflectance in the MIR domain. This surface type is restricted to narrow zones on the ice sheet located next to the ice sheet margin and corresponds to debris-covered ice.

Verification of glacial geomorphological mapping, accuracy assessment

Pixel-based accuracy assessment

The glacier-hydrological map of Thomsen and Winding (1988), chosen to assess the accuracy of the classification, covers the northeastern part of the area (green polygon in Fig. 3). This 1 : 75 000 scale photogrammetric map is based on aerial photographs acquired on July 10, 1985. To facilitate the pixel-based comparison, the hardcopy map was scanned and re-projected to the same coordinate system as the Landsat image. The main purpose of the map was to depict glacial geomorphological features, therefore features on the land area were generalized. For example, small patches of exposed rock and many small lakes and streams, clearly visible on the Landsat image, are absent in the published map. Trimline zones were only marked along the ice sheet and outlet glaciers, but not in the large valleys, formerly occupied by glacier tongues or outlet glaciers. We used information from historical photographs (Wedick, 1968; Mikkelsen and Ingerslev, 2002) to estimate the extent of the trimline zone in these valleys (Fig. 7).

Since the map did not distinguish between different vegetation types and water conditions, we could only perform a rigorous accuracy assessment on the classification map with the combined classes (Fig. 4: lower panel). Additionally, we combined the (d) debris-covered ice and (e) trimline zone classes because (1) the boundary between these classes was not mapped accurately in large areas (marked as “inferred” on the map), and (2) changes in the ice sheet boundary between the time of acquisition of the aerial photographs and the Landsat imagery were most significant in areas where a wide band of debris covered the ice sheet margin in the published map.

Four hundred random points, 100 points per class, were chosen from the map using a stratified random selection procedure (Richards and Jia, 1999: p. 267). To reduce the effect of the generalized nature of the map and the coregistration errors, a window size of 5 by 5 was used in the selection of the points, and a threshold value of 18 (out of 25) was stipulated before a point was eligible for selection. Points retained for the accuracy assessment were distributed among the classes as follows: ice, 98 points; water, 94 points; land, 81 points and trimline zone and debris on ice, 56 points.

Table II shows the error matrix from the accuracy assessment, the accuracy totals and the Kappa-statistics. The overall classification accuracy of 93.3% with an overall Kappa statistics of 0.9105 is very good. Discrepancies between the classified pixels and the reference data were potentially caused by actual retreat of the ice sheet margin (ice/snow/firn class replaced by trimline zone or debris on ice class) or by the presence of large barren (non-vegetated) areas over the ice-free land, mostly around lakes or along streams (land pixels which were classified as trimline zones or debris on ice).

Checking completeness and accuracy of trimline zone mapped from Landsat imagery

In addition to the pixel-based accuracy assessment using the Thomsen and Winding (1988) map, we also evaluated the accuracy and completeness of trimline zone (classes M and N in the upper panel and class e in the lower panel of Fig. 4, also depicted as the dark grey zone in the left panel of Fig. 7) around the ice sheet over a larger area. The ground truth map, shown in the right panel of Figure 7, is compiled from the map of Thomsen and Winding (1988) in the north and a map from Weidick (1992) in the south. Weidick’s map is based on the Geological Survey of Greenland topographic map sheet (1 : 250 000), which was drawn on the basis of point determinations and aerial photographs, The trimline zone obtained by the classification of Landsat imagery (dark grey) agrees very well with the trimline zone mapped by traditional methods. Errors include the exclusion of trimline zone pixels (error of omission) on steep slopes that cast shadows, causing trimline zone pixels to be classified as water (Fig. 7: circled region located south of Sermeq Avannarleq). Inclusion of non-trimline zone pixels into the trimline zone class (error of commission) is observed in large valleys (i.e. Pâkitsoq), where some glaciofluvial sediments are classified as trimline zone. Some bare patches along rivers or around lakes or in upland areas are also classified as trimline zone.

Figure 8 shows the location of the trimline, as measured by a handheld GPS receiver during the 2003 field campaign. The planimetric accuracy of the GPS measurements was about 10 m RMSE. The average distance between the trimeline locations derived from the Landsat imagery and the corresponding GPS measurements is less than the Landsat pixel size (28.5 m), demonstrating the geometric and thematic accuracy of the classified Landsat image. Red dots in Figure 8 show the boundary between regions with some vegetation and with absolutely no vegetation within the trimzone. These measurements and the fact that the boundary between the non-vegetated (class m, orange) and somewhat vegetated trimzones (class N, yellow) is nearly parallel with the inferred boundary of the paleo ice sheet, suggest that classes M and N mark trimline zone surfaces which have been exposed for different lengths of time.

Pansharpening of Landsat imagery for glacial geomorphological mapping

The detailed map of Thomsen and Winding (1988) also provides physiographic information, including features related to surface hydrology of the ice sheet. Mapped features include supraglacial lakes and streams, moulins, crevasses, as well as lineaments influencing the drainage pattern. Part of this map, covering the area north of Pâkitsoq (indicated by a yellow box in Fig. 3), is shown in Figure 9, together with a Landsat ETM+ based satellite image-map and the supervised classification map. The satellite image-map is a pan-sharpened false color composite of bands 2, 3 and 4. A special image fusion technique, the Bovery transform (Gupta, 2003: p. 268), was used to merge three multispectral bands with the higher spatial resolution panchromatic band. The fused image has the spatial resolution of the panchromatic band, while preserving the spectral information of the multispectral bands. The satellite image map is a very rich source of geometric information. It shows many details of the physiography of the ice sheet surface such as crevasse fields, supraglacial lakes, and flow traces that can be used to infer drainage patterns. On the adjacent land, the image shows glacial geomorphological features. This information can be combined with the classification map which distinguishes surfaces based on composition or nature of the land cover (e.g. debris-covered ice, fresh snow, etc.). Together, this provides a more comprehensive picture of the landscape than can be obtained from the pan-sharpened image or the classification map alone.

Table II

**Class accuracy assessment report**

Discussion

The primary motivation for the present study is the need to place recent changes observed on Greenland outlet glaciers in the broader context of glacier changes since the LIA maximum, i.e. 1850. Recognizing the limitations inherent to traditional field-based geological approaches, a multispectral Landsat ETM+ scene was analyzed to map trimline zones in the vicinity of Jakobshavn Isbrae, west Greenland. By utilizing information contained in the six non-thermal spectral bands, imaged surfaces can be classified according to their spectral reflectance characteristics. As demonstrated in Figure 4 (upper panel), this classification encompasses more than the distinction between bare rock and vegetation cover that defines the trimline. Indeed, different snow and ice facies, as well as water bodies with different turbidity are distinguished, and distribution of glaciofluvial and fluvial sediments and trimline subzones is mapped.

As our accuracy assessment confirms, we can reliably distinguish between vegetation-covered areas and trimline zones. It is well known that vascular vegetation has a very strong spectral contrast in the NIR. As several recent studies demonstrated (Bechtel et al., 2002; Van der Veen and Csatho, this issue) the shape of the reflectance spectra for lichens is very different from that of other surface types or landcovers as well. Notably, the NIR maximum, around 1 600 nm, is a unique feature of the lichen spectra. Since vascular vegetation is also characterized by distinct spectral signatures we believe that trimlines can be mapped reliably over different rock types using images that cover the broad spectral range up to 2 400 nm. As such, the techniques applied in this study present an advantage over trimline mapping from aerial photographs, which may be efficient where the bedrock consists primarily of gneiss, but which may not work very well in the case of dark, basaltic bedrock (Weidick, 1968).

In most cases, the trimline zone and debris-covered ice can be distinguished, but there may be problematic areas, depending on the distribution of surface material. The spectrum of debris-covered ice is a mixture of the spectrum of debris material and that of ice. Clearly, in places where the ice is fully covered with debris the boundary between debris-covered ice and ice-free morainal plains can not be detected by spectral analysis only. However, even relatively small patches of exposed ice reduce the reflectance in longer wavelengths.

In this study, supervised surface classification was used in which training samples of different cover types are selected by the operator. This is a feasible approach when considering small regions, as was done here, and it allowed assessing the usefulness of multispectral imagery for trimline mapping without the need for designing more automated techniques. For larger projects, however, automation of trimline mapping is not only possible but also desirable. Such automation can be achieved for example by using a multi-dimensional edge detector algorithm to find simultaneous changes in all or most spectral bands. This procedure would extract all boundaries, not only the transition between exposed rock and lichen and vegetation cover. Thus, to obtain the actual trimline, additional constraints are needed, possibly based on a geometric-topological model. For instance, defining the trimline zone as a narrow band bordering the ice margin, candidate boundaries for the trimline can be selected from the entire ensemble of boundaries. It should be noted, however, that in all likelihood the final selection will require some input from the operator.

Mapping trimlines provides information on the spatial extent of ice margin retreat. While horizontal changes can be determined reliably from Landsat imagery, obtaining a fully three-dimensional reconstruction requires additional information. Accurate determination of the elevations of trimlines is essential for estimating loss in ice volume. Printed maps are available for the region but their accuracy is insufficient for use as reference DEM. A better approach would be to use stereo pairs of aerial photographs to extract ground control points and to measure a DEM of the study region. Using these ground control points and the DEM, the Landsat image and the aerial photographs can be precisely co-registered and the elevation of the trimline be directly measured from the stereo photographs. An alternative solution is to use ASTER imagery. The ASTER sensor simultaneously provides multispectral imagery and stereo coverage for measuring surface elevations, allowing the user to measure the location and elevation of the trimline from the same data set. Moreover, fusion of spectral information with surface elevation provided promising results for mapping debris-covered glaciers (Paul et al., 2004), a problem very similar to the detection of the boundary between the trimline zone and debris-covered ice. We are currently in the process of conducting a pilot study for mapping trimline zone around the Jakosbhavns Isfjord from a combined data set of Landsat imagery and aerial photographs and from ASTER imagery.

In addition to mapping different types of surfaces based on their spectral properties, Landsat imagery can also be used for morphologic mapping of the ice surface, in particular features such as flow lines, surface lakes and streams, and crevasse zones. By merging three multispectral bands with the higher spatial resolution panchromatic band, the spatial resolution increases from 30 m to 15 m for the pan-sharpened image. This dramatically improves the ability to identify morphological features than can be used to describe and constrain the spatial pattern of ice flow and drainage.

In summary, the pilot study presented here has demonstrated that glacier trimlines can be mapped with confidence from multispectral Landsat images. The key to success is the wide spectral range of Landsat, from 400 to 2 400 nm. The procedures are sufficiently expedient to allow for mapping trimlines around the entire perimeter of the Greenland Ice Sheet. Doing so would provide a benchmark for its maximum extent during the LIA. Comparison of this benchmark with the present-day ice sheet would provide more quantitative estimates of changes in ice volume over the past ~150 years, and hence of the contribution of the Greenland Ice Sheet to global sea level since the Industrial Revolution. This would allow the estimated current contribution to be placed in a broader context and provide better understanding of the behavior of the ice sheet and its interaction with climate.

Trimline Mapping from Multispectral Landsat ETM+ ImageryCartographie de la limite maximale des glaces par imagerie multispectrale Landsat ETM+

Abstract

Résumé

Introduction

Characteristics of lichens

Previous trimline studies

Feasibility study: Jakobshavn Isbrae