Topographic maps and remote-sensing images

Our study uses 39 topographic maps at scales of 1 :50000 and 1 : 100 000. These maps were derived from aerial photographs acquired between 1971 and 1980 by the Chinese military geodetic service. They were scanned, georeferenced and rectified with a kilometer grid. The map coordinates were expressed using the Albers conic equal-area projection. Glacier outlines, digitized using ArcMap software, formed the basic reference data for changes in the Koshi River glaciers. Glaciers in Nepal are documented in the Nepal glacier inventory (Reference Mool, Bajracharya and JoshiMool and others, 2001), derived from topographic maps of India and Nepal between the 1960s and 1980s. However, the period covered by the topographic maps is almost 20 years, rather than the short period of a few years that would allow glacier reference data of nearly the same time to be produced; in addition, the topographic maps information was missed. Our study also uses four Landsat MSS images acquired in 1976 and 1977 which cover a period of <1 month (Table 1), to generate a reference inventory of glaciers in Nepal that fall within the Koshi River basin by manual delineation.

Table 1. Overview of the satellite image data sources applied

We also use six Landsat ETM+ scenes with 1 5m resolution acquired in 2000, and Landsat TM data with 30 m resolution acquired in 2009/10, to monitor glacier extent in 2000 and 2009/10. Unfortunately, it was difficult to select available images in 2009/10 due to the snow cover in Koshi River, especially on SSH. Therefore, our study uses several Landsat ETM+ scan-line corrector off (SLC OFF) data with 15 m resolution, acquired in 2009 (Table 1), as supplementary data to check the glacier extent covered by snow in 2009. This work will help reduce the glacier extent errors caused by snow. Landsat images are provided by the US Geological Survey (USGS; http://glovis.usgs.gov" and the Global Land Cover Facility. The image data were orthorectified or rectified using the Universal Transverse Mercator (UTM) coordinate system referenced to the World Geodetic System 1984 ellipsoidal elevation (WGS84). All images were subsequently reprojected in the Albers conic equal-area projection to be consistent with the topographic maps described previously.

Currently, different criteria can be used to extract glacier boundaries from satellite images, such as those based on spectral ratio (band4/band5, band3/band5, normalized-difference snow index (NDSI)), supervised classification and unsupervised classification (Reference Li, Sun and ZengLi and others, 1998; Reference PaulPaul and others, 2000, Reference Paul, Kääb, Maisch, Kellenberger and Haeberli2002; Reference Zhang, Su, T. Li and ZhongZhang and others, 2001). Furthermore, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery combined with digital elevation model (DEM) and thermal infrared ASTER imagery can be used to delineate glacier debris-cover extent semi-automatically (Reference Paul, Huggel and KääbPaul and others, 2004; Reference Shukla, Arora and GuptaShukla and others, 2010). Although the delineation methods were applied as accurately as possible, it is still difficult to detect glacier outlines correctly due to debris cover and snow cover. Consequently, manual digitizing of glacier extents was used for the multispectral images (MSS in 1976, Landsat ETM in 2000 and Landsat5 TM/Landsat7 ETM) in 2009/10.

Meteorological data

Meteorological data in China are obtained from the China Meteorological Data Sharing Service System http://cdc.cma.gov.cn/home.do Two weather stations (Tingri and Nyalam) were used to analyze climate changes during 1972–2010. Meteorological data in Nepal are provided by the International Centre for Integrated Mountain Development (ICIMOD; Reference RajchalRajchal, 2006). They span the period 1988–2008; some years have data missing.

Precision evaluation

Our study uses different spatial and temporal resolution data sources (Table 1). Hence, it is necessary to assess the uncertainty of the results to ascertain the accuracy. Numerous factors affect the certainty of glacier delineation in mountain regions (e.g. clouds and snow cover; the sensor’s spatial resolution; spectral characteristics; mountain shadows; mixed-pixel problems; image-processing problems; and disagreements between manually and automatically derived outlines). Glacier mapping uncertainty depends mainly on the spatial and radiometric resolution of the image sources and the conditions at the time of acquisition (Reference Bhambri, Bolch, Kawishwar, Dobhal, Srivastrava and PratapBhambri and others, 2013); less than half-pixel accuracy is usually achieved (Reference BolchBolch and others, 2010).

The analyst’s misinterpretation of glacier imagery also leads to inaccuracy, though guidelines were developed within the framework of the Global Land Ice Measurements from Space (GLIMS) initiative (Reference Raup and KhalsaRaup and Khalsa, 2007). For example, lateral moraine and end moraine should not be included as part of the glacier, while debris-covered ice should be included, but neither of these conventions is followed (Reference Nie, Zhang, Liu and ZhangNie and others, 2010). However, assessment of the uncertainty is difficult because images are not usually acquired at the same time or with the same spatial resolution (Paul and others, 2013).

Previous experience suggests that debris cover and spatial resolution are two key causes of uncertainty (Reference Williams, Hall and ChienWilliams and others, 1997; Reference Silverio and JaquetSilverio and Jaquet, 2005; Reference PaulPaul and others, 2013). The uncertainty associated with resolution concentrates on the edge of the glacier and the region of intersection between rock and glacier, so it is evaluated by analyzing the glacier’s edge pixels (Reference WangWang, 2013). We used a buffer area with half the sensor resolution to calculate the uncertainty in the software package ArcMap. The results show that the higher the sensor resolutions, the smaller the uncertainty. The glacier outlines derived from Landsat ETM+ in 2000 thus have a smaller uncertainty than those derived from Landsat MSS in 1976 and Landsat TM in 2009. Furthermore, the larger the glacier, the smaller the uncertainty. Glaciers less than 1 km² typically have an uncertainty larger than 8.2%; the maximum uncertainty is 33%, derived from Landsat MSS. The uncertainty associated with our analysis is shown in Table 2.

Table 2. Uncertainty of glacier area according to glacier area classification

Standard uncertainty propagation was used to calculate the overall uncertainty during 1976–2009:

(1)

where S is the overall uncertainty in this study and S _i is the uncertainty component at different times

The uncertainties estimated for the time periods 1976-2000, 2000-09 and 1976-2009 are 5.1%, 3.1% and 5.6%, respectively. The overall uncertainty in 1976, 2000 and 2009 is 5.8%.

Inventory of glaciers in the Koshi River catchment

Our analysis shows that there were 2061 glaciers with a total area of 3225 ±90.3 km² in the Koshi River catchment area in 2009 (Fig. 1). About 1282 of these glaciers (55.1% of the total number) are <0.5 km² in area, and these small glaciers make up 6.2% of the total glacier area. About 442 glaciers (20.1% of the total number) are in the area range 1-5 km², and these make up 23.3% of the total glacier area. There are 29 glaciers (1.6% of the total number) in the >20km² area class and these make up 40.6% of the total glacier area. Thus, although small glaciers are more numerous than large glaciers, a high proportion of the glacier area is constituted by the few large glaciers.

Figure 2 depicts the distribution of glacier population and area in different area-size classes normalized to the total number and area in the study region. The glacier number and area are larger on NSH than on SSH. There are 1290 glaciers on NSH, accounting for 61.7% of the total number of glaciers in the Koshi River catchment, with an area of 1961.0 ±54.90 km², accounting for 60.4% of the total glacier area; and 771 glaciers on SSH, accounting for 38.3% of the total number of glaciers, with an area of 1264.1 ± 35.4 km², accounting for 39.6% of the total glacier area. The mean glacier sizes on NSH and SSH are 1.52 and 1.64 km², respectively.

Fig. 2. The normalized frequency distribution of glaciers and the normalized distribution of their areas in the seven glacier area-size classes on SSH and NSH.

Glacier area change

Comparison of the total area of all glaciers in the Koshi River basin in 1976 with their total area in 2000 and 2009 (Table 3) reveals a reduction by -775.4 ±224.0 km², or 19.4 ± 5.6% (0.59 ± 0.17% a–¹). The number of glaciers also decreased, by 145, during the period 1976-2009. The mean size of glaciers in the study area decreased from 1.81 km² to 1.56 km². During the sub-period 1976-2000, total glacier area reduced by about 641.7 ±204.0 km², or 16.0±5.1%, with a mean rate of glacier area loss of 0.67 ±0.21% a–¹. During 2000-09, the glacier area decreased about 133.7 ±107.5 km², or 4.0 ±3.2%, with a mean rate of 0.44 ±0.36% a–¹. This suggests that the rate of glacier shrinkage has slowed.

Table 3. Glacier changes in Koshi during 1976–2009

Comparison of glacier changes on SSH and NSH shows that their total glacier areas shrank by about 321.3 ±88.8 km² and 454.2 ± 135.3 km², respectively, or 20.3 ± 5 . 6% (0.62±0.17%a–¹ ) and 18.8±5.6% (0.57±0.17% a–¹), respectively, during 1976-2009. During 1976-2000, glaciers on SSH and NSH shrank about 238.1 ±80.9 km² and 403.6 ± 123.2 km², or 15.0 ±5.1% (0.62 ±0.21% a–¹) and 16.7 ±5.1% (0.70 ±0.21% a^–1), respectively. There are some differences in glacier change rate during 2000-09. Over this period, glacier area on SSH decreased by 6.2 ±3.2% (0.68±0.36%a–¹) while glacier area on NSH decreased by 2.5 ±3.2% (0.27±0.36% a–¹). Consequently, glaciers on SSH retreated faster than those on NSH. The number of glaciers on SSH increased by about 8, but the number of glaciers on NSH decreased by 27.

Since glaciers with size <1 km² have a smaller total area and larger uncertainty in area (Table 2), we performed the glacier-area change analysis again, with their changes excluded. We found that the maximum contribution to glacier area loss comes from those glaciers within the area range 1-5 km², which contributed 32.0% of total area loss between 1976 and 2009; the second largest contribution comes from glaciers larger than 20 km² which contributed 18.6%. The glacier area loss between 2000 and 2009 accounts for 52% and 24% of the total area loss during the period 1976-2009 on SSH and NSH, respectively. Glaciers larger than 5 km² account for 27.0% of the area loss on SSH. We observed that most debris-covered glacier tongues larger than 10 km² on SSH do not retreat as much as other glaciers because the debris keeps them stagnant or limits their recession. Nevertheless, glaciers on SSH retreat at a higher rate (0.68 ±0.36% a–¹) from 2000 to 2009 because mass loss is focused in the clean ice or middle parts of glacier ablation zones.

One glacier (No. G086210E28316N), with an area of 7.44 km² in 1976, was divided into two glaciers in 2000, which had areas of 3.26 and 2.80 km² respectively, a total area loss of 18.5%. Both glaciers continued to shrink, yielding a decrease in area of 0.1 km², or 1.7%, between 2000 and 2009. The glacier with an area of 3.26 km² was almost instable. Another glacier (G086264E28309N) which was close to glacier No. G086210E28316N also underwent shrinkage. While both glaciers G086210E28316N and G086264E28309N retreated from 1976 to 2009, the glacier moraine lakes expanded where they were parts of both glaciers in 1976 (Fig. 3).

Fig. 3. Changes of example glaciers during 1976–2009 (background: Landsat ETM+ acquired in 2000, path 141, row 40). Black line: glacier outlines in 1976; red line: glacier outlines in 2000; blue line: glacier outlines in 2009.

Glacier changes in surrounding areas

Glacier retreat has been reported for many parts of the Hindu Kush-Himalaya against the background of worldwide warming (Reference Fujita and NuimuraFujita and Naimura, 2011; Reference BolchBolch and others, 2012; Reference Frey, Paul and StrozziFrey and others, 2012). However, these previous studies did not consider a whole basin nor did they offer a comparison between SSH and NSH. Our results establish that the loss of glacier area in the Koshi River catchment (0.59 ± 0.17% a–¹) is one of the fastest reported in the Himalaya during 1976-2009. For example, existing studies derived from 466 glaciers in Chenab, Parbati and Baspa basins, East Himalaya, show that glacier area reduced from 2077 km² in 1962 to 1628 km² in 2004, accounting for 21% of total glacier area (Reference KulkarniKulkarni and others, 2007). Glacier changes in West Himalaya were reported to be ~.4% shrinkage (0.33% a–¹) from 1976 to 2003 (Reference Ye, Yao, Kang, Chen and WangYe and others, 2006). Reference Mool, Bajracharya and JoshiMool and others (2001) found that the area of 153 glaciers in the Poiqu River basin decreased ~2.28 km² from 1988 to 2000, accounting for 5% (0.41%a–¹) of the total area. A total of 999 glaciers in the Pumqu basin with an area of 1461.84 km² decreased 10% (0.30% a–¹) in area and 9% and others’ and Jin and others’ results reflected only part of the glacier change in the Koshi River basin. Similar work was done in the Qomolangma region, where the glacier recession was 0.30% a–¹ during 1974-2008, with slower glacier retreat since 2000 (Reference Ye, Zhong, Kang, Stein, Wei and LiuYe and others, 2009). Reference Nie, Zhang, Liu and ZhangNie and others (2010) surveyed glacier changes in the Qomolangma nature reserve where glacier area decreased 502 ±0.04 km², or 15.6%, during 1976-2006. However, this result cannot be compared with our glacier change estimates for the Koshi basin, since these authors excluded potentially debris-covered regions from their glacier mapping.

Bedrock outcrop

Glacier recession usually occurs in the glacier tongue and glacier lateral where the end moraine and lateral moraine are formed. We also noted that glacier shrinkage occurs in a few parts of the accumulation region on SSH (Fig. 4). Due to steep terrain, conditions for ice development here are poor. These regions were sometimes snow-covered or had a small ice depth. However, the ice dies easily away and there is bedrock outcrop. This is part of the glacier accumulation region and bears the ice redistribution. In this work, we also consider the bedrock outcrop as glacier shrinkage. This is one reason why the glacier change is larger on SSH than on NSH from 2000 to 2009.

Fig. 4. An example of bedrock outcrop: (a) background: Landsat ETM+ acquired in 2000 (path 140, row 41); (b) background: Landsat TM acquired in 2009 (path 140, row 41). Both positions are bedrock outcrop covered by snow/ice in 2000.

Reasons for glacier variation

Glacier variation is a complex response to climate change mediated by a time lag which is affected by several factors (e.g. glacier size, glacier bed slope, glacier type). The quantitative relationship between glacier shrinkage and climate change is difficult to pin down. Tingri station, located on the highest-elevation site on NSH and close to the region where glaciers are distributed, provides a relevant climatic record for our study area. A warming trend in this station’s temperature record during 1972–2011 can be seen in Figure 5. Temperature has risen by 1.5°C in the past 39 years, with a sharp rise around 1997. Winter temperature has increased significantly (Reference Yang, Zhang, Qin, Kang and QinYang and others, 2011).

Fig. 5. Changes in mean air temperature and precipitation during 1972–2011 recorded at Tingri station.

The Tingri station precipitation record showed a slight increase during 1972–2011. There were dry years between 1980 and 1998 (Fig. 5). A separate record derived from ice cores located in the region indicates that precipitation decreased rapidly from the 1950s to the 1960s, then kept relatively steady, and decreased between 1980 and 1997 (Reference Ren, Qin, Kang, Hou, Pu and JingRen and others, 2004). Precipitation has increased again since 1997. Temperature increase generally increases glacier ablation rates, and reduced precipitation before 1997 led to less net ice accumulation in this period. The warming trend on SSH is consistent with NSH (Reference WangWang, 2013). Thus, we believe that glacier retreat during 1976–2000 can be attributed primarily to air-temperature increase and precipitation decrease. We further suggest that glacier retreat during 2000–09 was controlled by temperature increase, even though precipitation increased at the same time; the increased accumulation could not offset ice losses to surface melting driven by higher temperature.

In this study, we have presented a new glacier inventory for the whole Koshi River basin derived from data sources including multispectral satellite imagery. The new inventory indicates that there were 2061 glaciers with an area of 3225.1 ±90.3 km² in 2009. Comparison with serial glacier mapping derived from different data sources reveals that the glacier area in Koshi River basin shrank by 19.4±5.6% (0.59±0.17%a–¹) during 1976-2009. The number of glaciers has decreased from 2206 to 2061. Several glaciers (<1 km²) have vanished and a couple of glaciers were dismembered and divided into two or more individual glaciers during 1976-2009 (e.g G086210E28316N and G086264E28309N in Fig. 3). Glacier reduction was slightly faster on SSH (20.3 ±5.6%) than on NSH (18.8 ±5.6%) from 1976 to 2009, and especially from 2000 to 2009.