Introduction

Alpine glaciers, especially those in temperate zones, are regarded as one of the best natural indicators of climate change and they have generally receded during the 20th century (Reference Dyurgerov and MeierDyurgerov and Meier, 2000). Hence there is a need for detailed global monitoring of glaciers (Reference OerlemansOerlemans, 1994; Reference Haeberli, Frauenfelder, Hoelzle and MaischHaeberli and others, 1999; Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004). Glaciers on the Tibetan Plateau play an important role in the global climate system (Reference Zheng and ZhuZheng and Zhu, 2003). Researchers have focused on alpine glacier variations in high Asia over the last few decades and found a general reduction in glacier area in the region (e.g. Reference Mayewski and JeschkeMayewski and others, 1979; Reference MillerMiller, 1984). According to many studies (Reference Chen, Liu and JinChen and others, 1996; Reference SuSu and others, 1999; Reference Pu, Yao, Wang, Ding and ZhangPu and others, 2001; Reference Wang and LiuWang and Liu, 2001; Reference Jing, Ye, Jiao and YangJing and others, 2002; Reference Liu, Shen, Sun and LiLiu and others, 2002; Reference Lu, Yao, Liu, Ding and LiLu and others, 2002), glacial retreat varies spatially across the Tibetan Plateau. Glacier recession rates are lower inland and higher at the margin of the Tibetan Plateau (Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004).

Due to the large number and remoteness of most alpine glaciers, remote-sensing satellite techniques, including microwave data and optical imagery, have been frequently used in global-scale surveys. Landsat imagery, including the Landsat Multispectral Scanner (MSS; four spectral bands in the visible/near-infrared (VNIR) parts of the electromagnetic spectrum with 57 m resolution), the Landsat Thematic Mapper (TM; seven spectral bands from the visible to the thermal-infrared part of the spectrum with 28.5 m spatial resolution in VNIR) and the Landsat Enhanced Thematic Mapper Plus (ETM+; eight discrete bands with 14.25 m spatial resolution in the panchromatic band), has been one of the primary data sources for glaciological research (Reference Bindschadler, Dowdeswell, Hall and WintherBindschadler and others, 2001) because it provides glacier information in remote areas since the beginning of the series in 1972 (Reference Meier, Freden, Mercanti and WittenMeier, 1973). Imagery from the Advanced Space-borne Thermal Emission and Reflection Radiometer (ASTER), which employs 14 discrete bands with three bands in VNIR with 15 m resolution and a 1 5m resolution NIR along-track stereo band looking backward 27.68 from nadir (Reference KääbKääb, 2002), is also widely used for assessment of glacier dynamics in many programs (e.g. the Global Land Ice Measurements from Space (GLIMS) project (Reference BishopBishop and others, 2004)).

Naimona’nyi, the highest peak of the western Himalayan mountains, with an elevation of 7694 m, is located in the southwestern region of the Tibetan Plateau (30˚04’– 31˚16^’ N, 818–81˚47^’ E). Many researchers have focused on alpine glacier variations in the western Himalaya during the last few decades (e.g. Reference BishopBishop and others, 2004), yet little is known about glacier variations in the Naimona’nyi region. In this work, glacier variations in the Naimona’nyi region have been surveyed using a series of digital images (Landsat MSS in 1976, TM in 1990 and 1999, and ASTER in 2003; Table 1) and 1:50000 topographic maps produced from aerial photographs in 1974.

Table 1.

Digital satellite images used in the paper

Methodology

This paper compares the results of two alternative methods for studying glacier variations. The most popular methods of research on glacier variations by remote sensing focus on how to extract glacier information by band algebraic operation, digitization or classification using individual images from different times and then comparing the results from individual images. In this paper, we report glacier area growth and shrinkage in the Naimona’nyi region during the period 1976–2003 by comparing individual images. The weakness of this method is that it overlooks discrepancies that always exist among sequential images due to different seasons, resolutions and data sources. Thus, this paper develops a new method for studying glacier variations by means of Geographic Information System (GIS) and remote-sensing techniques. We integrate all classification results from individual images (1976, 1990, 1999 and 2003) by the algebraic operations in the Arc/Info grid module which enables us to determine mismatched information or unreasonable changes in glacier variations. The results obtained using both the ‘individual image’ and the ‘integrated’ method are presented here.

Multitemporal and multi-source digital satellite images have to be accurately orthorectified by a digital elevation model (DEM) before calculating changes based on pixels. Orthorectification eliminates the effects of perspective distortion and reduces the relief displacement on remotely sensed data. For the Naimona’nyi study region, we use the highest-available-resolution DEM data (1 : 50 000 scale, DEM5, cell size 25m). The horizontal accuracy of DEM5 with respect to the 1 : 50000 topographic map of the region is within 1.0 gridcell, i.e. 25 m.

Before orthorectification of the sequential satellite imagery, the height accuracy of the DEM was evaluated by comparing 331 elevation check points on the 1 : 50 000 topographic maps with the corresponding height values for the same locations in the DEM. We obtained an average height difference of 12.37 m, with a standard deviation of 18.52 m. Residual root-mean-square error (rmse) (Reference Hall, Bayr, Schöner, Bindschadler and ChienHall and others, 2003; Reference Stevens, Garbeil and Mouginis-MarkStevens and others, 2004; Reference KääbKääb, 2005) of the DEM with respect to the topographic map is +20.93 m. Maximum height deviations were –180.9 and +27.8m. The accuracy of orthorectification is within one image pixel. All image data used were converted from disparate sources to a common format defined in Arc/Info grid with Transverse Mercator projection and Krasovsky 1940 spheroid. Co-registration for all orthoimages is based on the 1:50000 topographic map (which is used as the common ‘base’), and all the co-registration errors are within one image pixel (Table 2).

Table 2.

Rmse of verification points in each co-registered image compared with the 1:50000 topographic map

Glaciers are mapped by unsupervised classification using band arithmetic, TM4–TM5 for Landsat imagery and TM3N– TM4 for ASTER imagery. This is a simple, robust and accurate method for glacier classification because of the very low reflectance of ice and snow in the near and middle infrared (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). All results are resampled to the same 30 m gridcell resolution, and some snow-covered non-glacier areas are removed by filter methods and manual editing. Classification accuracy of individual imagery series since 1976 was analyzed by Kappa technology (Reference CohenCohen, 1960; Reference CongaltonCongalton, 1991) using 260 random points. The KHAT accuracy is 91.2% from the results of the Kappa analysis.

A classification scheme was used in which non-glacier areas were assigned a single-digit value of 5, and glacier areas were assigned a single-digit value of 7. This classification scheme facilitates the algebraic operations among classification results from images in the Arc/Info grid module. The one-digit value (i.e. 5, 7) from each gridcell was integrated over the four classification results from images (from 1976, 1990, 1999 and 2003) by map algebra to generate a map with a four-digit value in each gridcell that simulates glacier change during the period 1976–2003 (Fig. 1). This enabled us to track glacier variations during the corresponding period both on maps and in tables. Instances in which gridcells indicate changes that occur more rapidly than is considered possible (e.g. 5757, 7575, etc.) are considered mismatches (e.g. 5757 stands for gridcell variation in ‘non-glacier area, glacier area, non-glacier area, glacier area’ in the 1976, 1990, 1999 and 2003 images, respectively). ‘Noise’, which is caused by mismatched information or unreasonable changes, such as different data sources or seasonal differences in snow cover in non-glacier areas, is eliminated by reclassification using the remap table (Table 3). All the results are recoded by the Reclass function using remap tables in the Arc/Info grid module. The four-digit integrated value from each gridcell (hereafter called the ‘integrated unit’) is used to identify advancing glacier areas, where non-glacier areas become glacier areas, and retreating areas, where glacier areas become non-glacier areas during different periods. These regions can be identified by both the integrated map (Fig. 1) and tabulated information (Table 4).

Fig. 1.

Glacier area variations in the Naimona’nyi region, 1976–2003.

Table 3.

Remap table of reclassification for reconstructing grid units of glacier variation process

Table 4.

Information table of the integrated grid, 1976–2003, by value of gridcells

Results

From the classification results from individual images, the area of glaciers in the Naimona’nyi region was 87.04 km² in 1976 and decreased to 79.39 km² in 2003 (Table 5). This shows an obvious decrease in glacier area, and the rate of change varies during different periods. Recession was 2.59 km² during 1976–90 (or 0.19 km² a⁻¹ on average), 0.80 km² during 1990–99 (or 0.09 km² a⁻¹ on average) and 4.27 km² during 1999–2003 (or 1.07 km² a⁻¹ on average). From the individual images the total area decrease between 1976 and 2003 is 7.66 km² or 8.8% (Table 5).

Table 5.

Glacier area change, 1976–2003, by individual image classification results

By the integrated method that we use in this paper, the glacier area was 84.41 km² in 1976 and 77.29 km² in 2003 (Table 6), a change of 7.12 km² or 8.4%. Accelerated recession of glacier area is clearly shown in Table 6. There is a recession of 2.37 km² during 1976–90 (or 0.17 km² a⁻¹ on average), 1.67 km² during 1990–99 (or 0.19 km² a⁻¹ on average) and 3.08 km² during 1999–2003 (or 0.77 km² a⁻¹ on average) (Fig. 1; Table 6). This indicates that glacier retreat in the western Himalaya is dramatic, and has accelerated in recent years.

Table 6.

Glacier area change, 1976–2003, after integration and reconstruction of basic units by remap table

Glaciers in the Naimona’nyi region have both advanced and retreated during the period 1976–2003. However, the area of glacier recession is much larger than the area of glacier advance (Table 3), and the advancing area decreases through time, while the retreating area increases.

Most of the area of glacier retreat occurred at the termini of glaciers in the southeast of the region (Fig. 1), while most of the area of advance occurred at the termini of glaciers in the northwest of the region. Glaciers in the southeast are retreating faster than those in the northwest.

Discussion

This paper explores an alternative method for studying glacier variations. Comparison of the results of glacier area from individual images (Table 5) with those from the integrated method indicates less variation in glacier area and a more obvious trend of accelerated retreat in the integrated and recoded grid data (Table 6) in the Naimona’nyi region. We could detect discrepancies in glacier extraction among the sequential images both graphically (Fig. 1) and numerically (Table 4). Since most of the ‘noise’ or mismatches could be detected and eliminated in the temporal spatial integrated map data by reconstruction using remap tables in the Arc/Info grid module, this provides a better quantification of glacier variations. Additionally, the integration method enables us to determine different variation characteristics during different periods over space in the research region. We were able to identify some misclassified areas (e.g. high-elevation rock ridges surrounded by glaciers that were misclassified as advancing glaciers due to snow cover; Fig. 1), something that is hard to determine by classification results from individual images only. These areas could be identified in the integrated map and can be verified in the near future using glaciological information or detailed field surveys.

Our research shows that the recession of glaciers in the Naimona’nyi region is more extensive and faster in recent years than before. The decrease of 8.4% of the total area during 1976–2003 is larger than the average 7% retreat of glaciers in high Asia since the 1960s (Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004). Alpine glacier recession is occurring in many places in Asia (e.g. 4.8% loss of the total area of 889 km² during 1969– 2002 in the Geladandong region (Ye and others, in press; 13.8% recession in the Ürümqi river drainage area during 1964–92 (Reference Chen, Liu and JinChen and others, 1996); 17% loss in the A’nyêmaqên mountains of the Yellow River source during 1966–2000 (Reference Yang, Ding, Liu, Lu and ChenYang and others, 2003); and 10.3% decrease in the western Qilian Shan during 1956–90 (Reference Liu, Shen, Sun and LiLiu and others, 2002)). The recession of glaciers in the western Himalaya has accelerated, which coincides with glacier variation trends in high Asia (Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004). The retreat is due to the negative glacial mass balance and is affected by rising temperature and decreasing precipitation over the Tibetan Plateau (Reference Yao, Wang, Liu, Pu, Shen and LuYao and others, 2004).

Glacier variations in the Naimona’nyi region during 1976–2003 also show spatial differences. Retreat usually occurred in the southeast, while advance always occurred in the northwest. However, the cause for advance or retreat of some glaciers in the region has not been identified because of the paucity of field-survey and regional meteorological data. Furthermore, this work only studied variations of glacier area. It does not include height variations which might show downwasting. In the future, we will extract and verify DEMs using various ASTER (bands 3N and 3B) data sources, and calculate vertical differences from DEM5 to study the volume variation of glaciers in the region. This will enable us to generate a more comprehensive view of glacier variations in the Naimona’nyi region, which will require detailed study in the near future.

Conclusions

The integrated unit with a four-digit value over the four images (from 1976, 1990, 1999 and 2003) in the Naimona’nyi region enables us to determine mismatched information or unreasonable changes in glacier variations. By integration of all sequential spatial data through time based on fundamental grid units, most of the mismatches could be detected and eliminated in the integrated data using remap tables in Arc/ Info grid. The comparison of glacier area from 1976 to 2003 (Tables 5 and 6) demonstrates that the integrated method did well in reducing discrepancies from various data sources, different seasons and different resolutions. The integrated method achieved more accurate results in glacier variations by applying the Arc/Info manipulations of the original classification results from individual images. This paper provides an alternative method of synthesized research on glacier variation using multi-source and multitemporal data.

Our results show that areas of both glacier retreat and advance exist in the Naimona’nyi region during the period 1976–2003; however, retreat dominates, and increases through time. Glaciers in the southeast are retreating more than those in the northwest. The 8.4% decrease of glacier area in the Naimona’nyi region is dramatic compared with other regions and the average glacier recession in high Asia since the 1960s.