Glacier changes on the Tibetan Plateau derived from Landsat imagery: mid-1970s – 2000–13

Abstract

Glacier area changes on the Tibetan Plateau were studied in different drainage basins based on Landsat satellite images from three epochs: 263 in the mid-1970s, 150 in 1999–2002 and 148 in 2013/14. Three mosaics (M1976, M2001 and M2013) with minimal cloud and snow cover were constructed, and the uncertainty due to each epoch having a finite span was accounted for. Glacier outlines (TPG1976, TPG2001 and TPG2013) were digitized manually with guidance from the SRTM DEM v4.1 and Google Earth imagery. To achieve complete multi-temporal coverage in a reasonable time, only debris-free ice was delineated. Area mapping uncertainty was evaluated at three study sites, Mount Qomolangma (Everest), Mount Naimona'Nyi, Mount Geladandong, where the largest differences between present and earlier measurements were within ~±4%. Area differences with previous inventories ranged from −19.6% (TPG1976 minus the first Chinese Glacier Inventory) to −3.6% and −1.1% (TPG2013 and TPG2001, respectively minus the second Chinese Glacier Inventory), while the difference TPG2001 minus the GAMDAM Glacier Inventory was +10.4%. Glacier area on the plateau decreased from 44 366 ± 2827 km² (1.7% of the study area) in the 1970s to 42 210 ± 1621 km² in 2001 and 41 137 ± 1616 km² in 2013. Shrinkage was faster in external drainage basins of the southeast than in the interior basins of the northwest, from a maximum of −0.43% a⁻¹ (−1.60% a⁻¹ during 1994–2013) in the Mekong catchment down to a minimum of −0.12% a⁻¹ in the Tarim interior drainage.

1. INTRODUCTION

Nearly all mountain glaciers have been melting rapidly in recent decades (Gardner and others, 2013; Zemp and others, 2015) and this has been documented, especially in High Mountain Asia (HMA hereafter), in surveys of the whole (Li and others, 2008; Cogley, 2016) or large parts (Bolch and others, 2012; Yao and others, 2012) of the region as well as in basin-scale and single-glacier studies (Fujita and Nuimura, 2011; Wang and others, 2013; Yang and others, 2013). As reported in these earlier studies, most HMA glaciers showed obvious shrinkage (area loss) and thinning (implying mass loss) except in the west, for example in the Karakoram (Gardelle and others, 2012) and western Kunlun (Neckel and others, 2014; Kääb and others, 2015; Ke and others, 2015). However, due to strong spatiotemporal differences of climate variability and topographic effects, the magnitudes of glacier shrinkage have varied in space and time. Thus it is necessary to quantify, in a unified scheme of change detection, not only the rates of shrinkage but also their spatial variations between different sub-areas (e.g. catchments) of the Tibetan Plateau (TP hereafter).

The TP, defined by outlines from Zhang and others (2014), can be divided into two categories of hydrological unit, the outlines and basin codes of which we take from Shi and others (2008). One is the interior-drainage basins (with no outlet to the ocean; IB hereafter), including the Central Asia IB, basin code 5Y and the TP IB, code 5Z (Shi and others, 2008). The Central Asia IB includes the internal drainages of Hexi (5Y4), Tarim (5Y6) and Qaidam (5Y5). The TP IB includes six sub-areas (5Z1–5Z6; Fig. 1). The other category is external river drainage catchments (EC hereafter), including those of the Yellow River EC (hereafter 5J), the Yangtze River EC (5K), the Salween River EC (5N), the Mekong River EC (5L), the Ganges River EC (5O) and the Indus River EC (5Q; Fig. 1). Generally, meltwater from glaciers in IB flows into inland depressions, while those in EC flow into the ocean and so contribute more directly to sea-level rise. Climatically, the Indian monsoon dominates in the southeastern TP in summer (e.g. in 5O, 5N, 5L), with some effects from the East Asian monsoon in the east (e.g. in 5Y4, 5J, 5K). Glacier changes in these regions are therefore controlled by the summer monsoon atmospheric circulation system (Yao and others, 2012). The westerlies are dominant in winter, especially in the western TP (e.g. 5Y6, 5Z4, 5Z3). Regions in the interior TP, such as 5Z1, 5Z2, 5Z5, 5Z6 and 5Y5 are less influenced both by the summer monsoon and the westerlies (Yao and others, 2012).

Fig. 1. Landsat data used for glacier outlines in this study. The Central Asia IB includes internal drainages of Hexi (5Y4), Qaidam (5Y5), and Tarim (5Y6). The TP IB includes six sub-areas (i.e. 5Z1–5Z6). The EC at the marginal TP include the Yellow River EC (5J), the Yangtze River EC (5K), the Salween River EC (5N), the Mekong River EC (5L), the Ganges River EC (5O) and the Indus River EC (5Q). Glaciers are in black.

In this study, we present debris-free glacier area changes and their spatial heterogeneity at the plateau scale from the 1970s to 2013 based on three mosaic Landsat images, i.e. in the 1970s (TPG1970s), TPG2001 and TPG2013.

Another purpose of the study is to enrich the available collection of glacier inventories in HMA, by presenting a detailed comparison of TPG datasets with previous glacier inventories. We compare the time series of TPG with two newly released glacier inventories: the Second Glacier Inventory of China, based on satellite images during 2004–11, 92% of which were obtained during 2006–10 (Version 1.0, CGI-2 hereafter) (Guo and others, 2014, 2015; Liu and others, 2015) and the GAMDAM (Glacier Area Mapping for Discharge from the Asian Mountains) glacier inventory for the period 1999–2003, GGI (Nuimura and others, 2015); and with the first glacier inventory of China, CGI-1 (Li and others, 2008; Shi and others, 2008) for the period 1956–84, which was based on more than 2000 1:50 000 and 1:100 000 topographic maps derived from more than 342 000 aerial photographs from the 1950s to the 1970s, and with reference to 200 Landsat MSS images in the 1980s and 1990s.

2. DATA AND METHODS

2.1. Data

We obtained the boundary data of the TP from Zhang and others (2014) and of the drainage basins from CGI-1 (Shi and others, 2008). The basin boundaries were also used by Guo and others (2015) for CGI-2. Marginal catchments, including all ECs (i.e 5O, 5Q, 5J, 5K, 5N, 5L) and two of the IB (5Y6 and 5Y4) were truncated at the TP outline. Together these catchments cover 2.57 × 10⁶ km² (1.46 × 10⁶ km² in the IB and 1.11 × 10⁶ km² in the EC), which is larger than one quarter of the area of China.

Time series of Landsat scenes were processed with Standard Terrain Correction (Level 1T), which provided systematic radiometric and geometric accuracy by incorporating ground control points while employing SRTM DEM for topographic accuracy. L1T satellite images were chosen since the 1970s (Table 1, Table S1) from the US Geological Survey at http://glovis.usgs.gov/. High-quality images with minimal cloud and snow coverage were selected for the delineation of glaciers. Both wintertime and summertime images were considered acceptable in the interest of minimizing seasonal snow and cloud.

Table 1. Satellite image list for TP glacier mapping

We treated 5% coverage by snow and clouds as the acceptable threshold, but were obliged to work with some images that had more than 20% snow or cloud cover. To minimize the effects of snow or cloud cover on glacierized areas, we made supplementary use of high-resolution (30 m spatial resolution and 4 d repetition cycle) images covering almost the whole TP from the Chinese satellites HJ-1A and HJ-1B, which were launched on 6 September 2008. Both carried as payload two 4-band CCD cameras with swath width 700 km (360 km per camera). All HJ-1A/1B data in 2012, 2013 and 2014 (65 scenes, Fig. S1 and Table S1) were from China Centre for Resources Satellite Data and Application (CRESDA; http://www.cresda.com/n16/n92006/n92066/n98627/index.html). Each scene was orthorectified with respect to the 30 m resolution DEM of the Shuttle Radar Topography Mission (SRTM) and Landsat images in the 2000s.

For the 1970s, 189 L1T scenes of Landsat MSS (60 m spatial resolution) were collected over the period 1972–79, covering 92% of the whole TP study area (Fig. 1; Table S1), including 116 scenes in 1976/77. However, for some areas without qualified images in the 1970s, 16 images of Landsat TM were acquired in later years (1981, 1986–89 and 1994). Of the TP study area, 6.5% was covered by the 1980s images (Table S1) and 1.3% by two 1994 images (Fig. 1).

A total of 150 Landsat7 ETM+ scenes (spatial resolution: 30 m) was selected for 1999–2002, including 61 scenes (41%) in 2001 and 47 scenes (31%) in 2000 (Table S1). For ~2013, 128 Landsat 8 Operational Land Imager (OLI) images were selected with 30 m spatial resolution for comparability with previous and current glacier inventories; ~20 images acquired in 2014 were used to complete the coverage of the TP. All satellite data were mosaicked to represent three periods, the 1970s, 2000s and 2013. The most frequent year in each period was defined as the reference year for the mosaic image: i.e. 1976, 2001 and 2013. The mosaic image in the 1970s (M1970s hereafter) was separated into three different sub-periods to enable us to calculate rates of change more accurately: M1970s-76, M1970s-88 and M1970s-94. For the same reason, the M2001 and M2013 mosaics were also spatially separated into the corresponding subsets M2001-76, M2001-88, M2001-94, and M2013-76, M2013-88 and M2013-94. All selected Landsat scenes are presented in Figure 1.

Large-scale topographic maps from aerial photographs taken in the 1960s and 1970s were used to validate glacier outlines in complex terrain, including 97 maps at 1:50 000 scale and 85 at 1:100 000 scale to obtain geometrical features at the three validation sites (Mount Qomolangma, Mount Naimona'Nyi and Mount Geladandong). The quality of these topographic maps was evaluated in the original studies of the validation sites by Ye and others (2006a, b, 2009), and most were found to be of good quality as basic reference data. However, there was some misinterpretation of snow cover as ice (Fig. 8).

2.2. Glacier outline delineation and quality control

Glacier outlines were digitized on-screen manually from the three image mosaics, relying on false-colour image composites (MSS: red, green and blue (RGB) represented by bands 321; TM/ETM+: RGB by bands 543, OLI: RGB by bands 654), which allowed us to distinguish ice and snow from cloud. We also referred to the SRTM DEM v4.1 and overlaid our glacier outlines on Google Earth imagery of different dates and times, which allowed us to minimize the effects of mountain shadows and seasonal snow cover.

Debris-free ice was distinguished from the debris and debris-covered ice by its higher reflectance. Debris-covered ice was not delineated in this study, as the positions of debris-covered glacier termini can be difficult to distinguish in the multispectral images. For brevity we refer below to ‘glacier area changes’, but it should be understood that this means ‘changes in the area of debris-free ice’. The implications are discussed in Section 4.

Adjacent glaciers were regarded as glacier complexes, i.e., as continuous unsubdivided ice bodies. This point is also discussed in Section 4.

The manual delineation was carried out by three operators over a period of 6 years. As the operator's experience and skill may affect the accuracy of delineation, the delineated glacier outlines were validated and edited by different operators at least three times by overlaying them onto Google Earth imagery, DEMs, topographic maps and corresponding satellite images by another two glacier researchers with field survey experience. Figure 2 presents an example of glacier outlines from TGP2001 in different versions edited by different operators, together with black outlines from CGI-2. It suggests that operator experience could affect the accuracy of glacier delineations. In this example the yellow line was accepted as the final mapping result.

Fig. 2. Examples of (a) glacier outlines from different operators and different methods; (b) Separation of debris-free and debris-covered ice. (The background image was from Landsat7-ETM+ on 30 October 2000).

We also made some comparisons with band-ratio (e.g. TM3/TM5; Paul and Kääb, 2005) results, as shown in Figures 2, 3. Compared with manually delineated glacier outlines, automated results can be inaccurate or misclassified in some situations, where in addition to glacier ice there are also glacial lakes, clouds or mountain shadows. Figure 3 presents four examples of glaciers in mountain shadows misclassified as non-glacier areas by the band-ratio method. These comparisons show that the image-specific selection of thresholds for the band-ratio method would have required time-consuming manual post-editing, leading us to prefer manual methods directly.

Fig. 3. Evaluation of TPG2001 mapping results with the band ratio ETM3/ETM5. (a) Background image on 30 October 2000 from Landsat7-ETM+. (b) Background image on 13 May 2015 from Google Earth).

In this study, for areas with mountain shadows and snow cover, we tried to verify different methods using data from different seasons. For example, as mountain shadows are generated at different solar elevation angles in different seasons, glaciers in mountain shadows and snow cover could be delineated in a TPG mosaic image in December, while the accuracy of delineation could be checked in images in a different season with fewer shadows (Figs. 4, 5). Similarly, for glaciers in deep shadow, Google Earth™ imagery from different dates was used as the reference for manual delineation (Figs. 3, 5).

Fig. 4. Seasonal snow cover in (a) September and mountain shadows in (b) December in 2014 (RGB:742).

Fig. 5. TPG2013 glacier outlines in deep mountain shadows. (a) Background image on 30 November 2014 from Landsat8-OLI. (b) Background image on 13 May 2015 from Google Earth).

Steep slopes or headwalls were also excluded from the glacier coverage data in the three TPG epochs. Topographic maps from the 1970s (e.g. Fig. 6) and all available satellite images (including Google Earth™ imagery and HJ-1A/1B satellite data) were used as base reference data. Areas that appeared in any of these sources to have the characteristics of exposed ground/basement/bedrock were manually delineated as non-glacier, and were also cross-checked with CGI-1 (Shi and others, 2008) and CGI-2 (Guo and others, 2015). Steep hanging glaciers were included in each TPG inventory if they were identifiable on images in all three epochs. However, some smaller glaciers that disappeared in later images were retained in the earlier TPG data.

Fig. 6. Steep slopes or headwalls were excluded by referring to topographic maps of the 1970s. (a) Excluding steep slope areas without glacier ice in 2013. (b) Steep slopes without glacier ice on topographic maps in 1974.

The accuracy of outlines was controlled during digitizing by keeping the cursor within one half pixel of the interpreted ice margin. All glacier areas were calculated on the WGS84 spheroid in an Albers equal-area map projection centred at 95°E, 30°N with standard parallels at 15°N and 65°N.

2.3. Area change analysis

Glacier area change was analysed in both physical (ΔA; km²) and normalised percentage (AP; %) units, and also as rates of change (AR; % a⁻¹):

(1) $$\Delta A = A_1 - A_0 $$

(2) $$AP = 100\,\Delta A/A_0 $$

(3) $$AR = AP/\Delta t.$$

Here Δt = t ₁ − t ₀ is the time span of the calculation in years; t ₀ is the earlier time in the comparison, for example, 1976, 1988 or 1994, and t ₁ is the later time, for example 2001 or 2013. Rates of change AR were calculated separately for the three subsets, M1970s-76, M1970s-88 and M1970s-94, of the M1970s mosaic (Tables 2–4).

Table 2. M1970s-76* glacier area changes in different drainages on the Tibetan Plateau (km²). Basins in bold font are first-order basins; those in roman font are second-order basins

* M1970s-76 was the 1970s satellite images in the 1970s mosaic. Glacier cover%, glacier areal percentage in the drainage. Glacier IB/EC %, glacier area contribution of all glacier areas in the overall IB or EC drainage.

Table 3. M1970s-88* glacier area changes in different drainages on the Tibetan Plateau (km²). Basins in bold font are first-order basins; those in roman font are second-order basins

* M1970s-88 was the 1980s satellite images in the 1970s mosaic.

Table 4. M1970s-94* glacier area changes in different drainages on the Tibetan Plateau (km²). Basins in bold font are first-order basins; those in roman font are second-order basins

* M1970s-94 was the 1994 satellite images in the 1970s mosaic.

2.4. Uncertainty assessment

As shown in Section 2.2, glacier data accuracy was limited by several key factors, including obscuration of glacier boundaries by seasonal snow remnants, mountain shadows and clouds, registration and co-registration of satellite images, manual digitization accuracy and differences between methods. As only debris-free ice was studied in this work, there are no errors arising from the delineation of debris-covered ice, although we discuss this point later.

The error in glacier area was assessed using glacier-outline buffers (Rivera and others, 2005; Granshaw and Fountain, 2006), as summarized in Table 5. The uncertainty of the outline position, chosen to be one half of an image pixel (Bolch and others, 2010), i.e. ±30 m for Landsat MSS and ±15 m for Landsat-7 TM/ETM+, Landsat8 OLI and HJ-1A/1B, was multiplied by the length of the margin. As HJ-1A/1B data were supplementary to Landsat8 OLI at the same pixel size for the same vector data, its glacier outline uncertainty was involved in the calculation in 2013 (Table 5). We assume there is no error in the digitized length, so the area uncertainty varies only with the pixel resolution of the images. We express the area uncertainties as percentages of total glacier area.

Table 5. Uncertainties of glacier coverage by percentage (U, %) and uncertainty buffer area (UA, km²) in different basins. Basins in bold font are first-order basins; those in roman font are second-order basins

Glacier area uncertainty was also evaluated at three regions studied earlier by different methods (Table S2): Mount Qomolangma (Ye and others, 2009), Mount Geladandong (Ye and others, 2006a) and Mount Naimona'Nyi (Ye and others, 2006b). Glacier area difference (AD) was calculated between different datasets at similar times; its calculation was similar to that of ΔA (Eqn (1)); however ΔA represents glacier area change between different epochs of the same dataset.

Cogley (2016) shows that uncertainty in the time span Δt and uncertainty in the area change ΔA contribute equally to the uncertainty in the shrinkage rate AR. To calculate the standard error of a date or date range, the number of days between the start and end of the range was summed and divided by 4 (an arbitrary choice), and converted back to decimal years (Cogley, 2016). For date ranges such as 1972–79, 7 years, this means an uncertainty of ±1.75 years, with 64% of the area covered by images from 1976/77 in M1970s-76. For date ranges such as 1981–89, 8 years, the uncertainty is ±2 years, with 61% of the area covered by images from 1987/88 in M1970s-88. For date ranges such as 1999–2002 and 2013/14, the uncertainty is ±1.0 years and ±0.5 years, respectively.

In our calculations, the uncertainty of any area change se(ΔA) was obtained by adding the area uncertainties in quadrature (Cogley, 2016):

(4) $${\rm se}\,(\Delta A{\rm )} = \sqrt {({\rm se}\,A_0 )^2 + ({\rm se}\,A_1 )^2} $$

and the uncertainty se(Δt) of any time span is

(5) $${\rm se}\,(\Delta t) = \sqrt {({\rm se}\,t_0 )^2 + ({\rm se}\,t_1 )^2} . $$

The uncertainty of the corresponding shrinkage rate se(AR) is

(6) $${\rm se}(AR) = \sqrt {(AR \times {\rm se}\,\Delta t)^2 + ({\rm se}\,\Delta A)^2} /\sqrt {({\rm se}\,\Delta t)^2}.$$

3. RESULTS

3.1. Accuracy of TPG data

Given the assumed width uncertainties of ±30 m in 1976 and ±15 m in 2001 and 2013, separated by different periods in the 1976 mosaic image (Fig. 1), the uncertainty of total glacierized area was ±6.4% in TPG1976, while the other two uncertainties were smaller, i.e. ±3.8% in TPG2001 and ±3.9% in TPG2013. The uncertainty in area changes was 7.4% for 1976–2001, 5.5% for 2001–13 and 7.5% for 1976–2013. The uncertainty in the time span is ±2.0 years for 1976–2001, ±1.0 years for 2001–13 and ±1.8 years for 1976–2013.

The uncertainty in the corresponding shrinkage rates is ±0.15% a⁻¹ for 1976–2001, ±0.24% a⁻¹ for 2001–13, and ±0.17% a⁻¹ for 1976–2013.

In addition, comparisons between TPG data and our previous studies show that the differences at all three sites between the series of TPG data and the local studies are within ±4% (Table S2). The largest difference (present study minus previous study), ~ −4%, was at Mount Naimona'Nyi. Part of the difference can be attributed to different acquisition dates (Table S2; Fig. 7a): 13 October 2001 for TPG2001 and 9 November 1999 in Ye and others (2006b). Two years of shrinkage at −0.77% a⁻¹ would reduce the difference by almost one quarter. Another large difference, −3.4%, was found in 1976 using the same Landsat2 MSS images between 1976TPG and the 1976 Mount Naimona'Nyi imagery (e.g. Fig. 7b). This difference was due to different operators and different methods, i.e. manual digitization in this work and an integrated method in Ye and others (2006b).

Fig. 7. Comparison between TPG data and previous studies. (a) TPG2001 in yellow line and 1999 Mount Naimona'Nyi glacier in red; background image on 9 November 1999 from Landsat7-ETM+. (b) TPG1976 in yellow and 1976 Mount Naimona'Nyi glacier in red; background image on 6 December 1976 from Landsat2 MSS.

3.2. Comparison among different glacier inventories

As the Randolph Glacier Inventory version 3.2 (Pfeffer and others, 2014) was taken from CGI-1 and version 5.0 was taken from CGI-2 over the study area (except for the southeastern TP, taken from GAMDAM), we did not compare it with our results. Calculations of glacier area difference (AD) and area difference rate (AR) were made as described in Sections 2.3 and 2.4.

For comparisons between TPGs and CGI-2, we reduced the total glacierized area given for each region in CGI-2 by subtracting the corresponding debris-covered area. For comparisons between TPGs and CGI-1, the total CGI-1 glacier area in each region was reduced by the percentage given for debris cover in CGI-1. According to the CGI-2, there are 41 649 single glaciers covering an area of 42 682 km² in our study area on the TP. A comparison between CGI-1 and CGI-2 has been conducted by Guo and others (2015). Therefore, we only present glacier area differences between the time series of TPG and the two CGI datasets (Table 4). We were unable to remove debris-covered areas from the GAMDAM inventory (GGI) (Nuimura and others, 2015), which includes debris cover but does not distinguish it from debris-free ice.

The area difference TPG1976 minus CGI-1 was −19.6%. At sub-regional scale, the greatest area difference was ~ −41% in sub-region 5Z3. Possible explanations for such large differences include differences between the dates of the Landsat MSS images selected for TPG1976 and the topographic maps and air photos that were the source for CGI-1, and mis-interpretation of glacier outlines affected by snow cover, clouds or mountain shadows. Such misinterpretation is possible in either inventory, but we note that Gardelle and others (2013) inventoried a study area in the southeastern TP and found a glacierized area in 2011 (1303 km²) that was much less than the CGI-1 area (2541 km²). The implied shrinkage rate of ~−1.4% a⁻¹ over ~35 years is implausibly large for this area and suggests that CGI-1 may have included too much snow or have involved misidentified glacier areas. Similarly, Bolch and others (2010) present evidence that the Chinese topographic maps misinterpret glacier outlines in their study area in Mount Nyainqentanghla.

Between TPG2001 and CGI-2, the area difference was ~ −1.1% (Table 4). Between TPG2013 and CGI-2, the area difference was −3.6%. The range of dates in CGI-2 is from 2006 to 2010, about 8 years later than TGP2001 and 5 years earlier than TPG2013, so glacier shrinkage presumably explains part of these differences. For example, from the re-compiled results of Cogley (2016), glaciers over HMA as a whole from 1960 to 2010 were changing at −0.40% a⁻¹, which suggests that CGI-2 should be ~3% smaller (instead of ~1% bigger) than TPG2001 and ~2% bigger than TGP2013 (instead of ~4% bigger). Considering the measurement uncertainties, the actual differences between TPG and CGI-2 thus seem to be consistent with each other and with an HMA-wide average shrinkage rate of −0.40% a⁻¹. Moreover, the discrepancies between these three independently compiled inventories represent good agreement by any external standard. For example they are smaller than the area uncertainties of ±7.7 to ±8.4% estimated by an independent method for the three regions of the Randolph Glacier Inventory that overlap the TP (Pfeffer and others, 2014).

In general, differences between CGI-1/CGI-2 and results from previous local sites were larger than differences from the TPG datasets. This is true especially for CGI-1, which was 25% larger than the local dataset in 1974 in Mount Naimona'Nyi region (Table S2). Figure 8 shows that the largest differences mainly originated from two sources. One is the debris-covered areas in the south-west part of the region, and the other is over digitized glacier areas, especially in CGI-1. The difference can be explained in terms of differences of interpretation of satellite images in the local study and topographic maps from aerial photographs in CGI-1. Figure 8 inset shows an example of mis-interpreted glacier area from topographic maps in the Mount Naimona'Nyi region.

Fig. 8. Glacier outlines in TPG/CGI-1/CGI-2 in Mount Naimona'Nyi region; background image on 6 December 1976 from Landsat2 MSS. Inset: an example of mis-interpreted glacier area in CGI-1 from a topographic map in 1974.

The GAMDAM Glacier Inventory (GGI) (Nuimura and others, 2015) covers all of HMA for the period 1999–2003. GGI has less glacierized area than TPG 2001 in all basins. The area difference for the whole TP, TPG2001 minus GGI, was 10.4% (Table 6). The largest difference was 23.2% in 5O region, the Ganges River EC. The next largest differences, 15.4% and 13.8%, were in the Yangtze River EC (5K) and the Mekong River EC (5L). For the same three regions, the differences between TPG and CGI-2 are more moderate; AP _T2001-CGI2 was +2.3% in 5K, +7.9% in 5L and −10.4% in 5O (Table 6). These contrasting differences are discussed at greater length in Section 4.

Table 6. Comparison between glacier inventories in China or HMA. Basins in bold font are first-order basins; those in roman font are second-order basins

A_CGI-1 and A_CGI-2, are clean glacier area that excludes debris-covered areas. AP_T1976-CGI1(%), is the normalized area difference between TPG1976 and CGI-1, which was calculated by: AP_T1976-CGI1(%) = (A_T1976 − A_CGI-1)/A_CGI-1 × 100%. AD_T2001-GGI (km²) is glacier area difference between TPG2001 and GGI, i.e. AD_2001−GGI(km²) = A_T2001 − A_GGI. ADI_CGI-2, is the area of debris-covered ice from CGI-2.

No._CGI1, No._CGI2, No._GGI is the number of glaciers in CGI-1, CGI-2 and GGI, respectively.

3.3. Glacier shrinkage on the TP

Glacier shrinkage on the TP was substantial but spatially variable between the 1970s and 2013 (Fig. 9). It was more rapid in the EC than in the interior IB. Glacier area and its change were analysed by basin (Fig. 10; Table 2). In the 1970s, glaciers covered ~1.7% of the TP area. The total glacier area on the TP decreased from 44 366 ± 2827 km² in the 1970s to 42 210 ± 1621 km² in ~2001 (area change was −2156 km²), and further to 41 137 ± 1616 km² in ~2013 (area change was −1073 km²). In total, glaciers on the TP have changed by −3229 ± 3256 km² during 1976–2013.

Fig. 9. Glacier area changes in Mount Geladandong area in different drainage basins; background image on 17 December 1976 from Landsat2 MSS.

Fig. 10. Glacier area changes in different basins. M1970s-76, M1970s-88 and M1970s-94 are subsets of the M1970s mosaic representing different time spans of satellite image coverage. M1970s-94 and M1970s-88 represent the satellite image coverage from 1994 and 1980s in the M1970s mosaic, respectively. M1970s-76, M1970s-88 and M1970s-94, stand for the characters in blue, brown and purple, which represent basin code and glacier area change rates in 1976–2013, 1988–2013 and 1994–2013, respectively.

Before presenting shrinkage rates in detail, we repeat the rate uncertainties given in Section 3.1: ±0.15% a⁻¹ for 1976–2001, ±0.24% a⁻¹ for 2001–13 and ±0.17% a⁻¹ for 1976–2013. Not all rates for different sub-periods (Section 3.3.1) and different basins (Section 3.3.2) are distinguishable reliably from zero, but all are most-probable estimates, and patterns that are consistent across time or space can be interpreted with greater confidence.

3.3.1. Glacier area changes for three different sub-periods in M1970s

Three different parts of the M1970s mosaic, as shown in Figure 1 and Tables 2–4, are considered here. Glaciers that were first imaged in the mid-1970s (M1970s-76) covered an area of 34 364 ± 2490 km², which represented 77.45% of the glacierized area. These glaciers shrank at −0.19% a⁻¹ up to ~2001 and at −0.20% a⁻¹ from ~2001 to ~2013.

Glaciers that were first imaged in the 1980s (M1970s-88) covered an area of 7462 ± 337 km², representing 16.8% of the glacierized area; 89% of the area was located in the EC (Fig. 1; Table 3). These glaciers changed by −0.41% a⁻¹ from ~1988 to ~2001, and by −0.27% a⁻¹ from 2001 to 2013. The rate for ~1988 to ~2013 was −0.34% a⁻¹.

Glaciers that were first imaged in 1994 (M1970s-94) covered an area of 2541 ± 337 km², representing 5.7% of the glacierized area. The entire area was located within the EC (Fig. 1; Table 4). These glaciers changed by −0.42% a⁻¹ from 1994 to 2013.

The main result of this comparison is that shrinkage rates clearly show an increase from the north-west to the south-east of the TP.

3.3.2. Glacier area changes in different basins

As shown more fully in Tables 2–4, glaciers in the EC tended to shrink more rapidly between the mid-1970s and ~2013, by −8.7% (−0.24% a⁻¹), than those in the IB (−6.2% or −0.17% a⁻¹). For the three sub-periods, shrinkage rates in the EC were −0.25% a⁻¹ (mid-1970s to ~2013), −0.32% a⁻¹ (~1988 to ~2013) and −0.42% a⁻¹ (1994 to ~2013). For the IB, the corresponding rates during the first two sub-periods were −0.16% a⁻¹ and −0.49% a⁻¹.

The most rapid shrinkage occurred in the southeast EC. Glaciers in the Salween basin (5N) shrank at −0.31% a⁻¹ (mid-1970s to ~2013), −0.36% a⁻¹ (~1988 to ~2013) and −1.04% a⁻¹ (1994 to ~2013), and those in the Mekong basin (5L) shrank at −0.31% a⁻¹, −0.43% a⁻¹ and −1.60% a⁻¹ during the same three sub-periods. This evidence for accelerating shrinkage seems clear, but it must be assessed with caution because the three sub-periods sample different glaciers and the sample sizes are comparatively small. Total glacierized areas in the mid-1970s (Table 6) were 1582 km² in 5N and 252 km² in 5L. Nevertheless the Ganges basin 5O, the EC with the largest glacierized area in the mid-1970s (13 965 km²), exhibits slower rates but a similar acceleration: −0.21% a⁻¹ (mid-1970s to ~2013), −0.28% a⁻¹ (~1988 to ~2013) and −0.35% a⁻¹ (1994 to ~2013). Shrinkage was also slower in the other ECs but with more limited evidence for acceleration. Rates for 5J (upper Yellow River) were −0.23% a⁻¹ for the mid-1970s to ~2001 and −0.26% a⁻¹ for ~2001 to ~2013, while the corresponding rates for 5K (upper Yangtze River) were −0.28% a⁻¹ and −0.39% a⁻¹ and for 5Q (upper Indus River) they were −0.32% a⁻¹ and −0.33% a⁻¹.

The Central Asia IB (5Y) accounts for 70% (17 100 km²) of the glacierized area of the IB (Table 2) and showed the slowest shrinkage: by −0.14% a⁻¹ for the mid-1970s to ~2001 and −0.13% a⁻¹ for ~2001 to ~2013. However, there is marked geographical variation within 5Y. In the northeast, 5Y4 (Hexi) and 5Y5 (Qaidam) considered together show more rapid shrinkage, especially for glaciers first imaged in the 1980s (Table 3), which shrank at −0.65% a⁻¹ before ~2001 but at only −0.35% a⁻¹ after ~2001. The slow rates for 5Y as a whole are dominated by the slower shrinkage in 5Y6 (Tarim) to the west, which accounts for more than four fifths of the total glacierized area of 5Y. (Sub-basin 5Y6 does not appear in Table 3 because all its glaciers were first imaged in the mid-1970s, i.e. they are part of M1970s-76.)

In the TP IB (5Z), in which all the glaciers are part of M1970s-76, the shrinkage rate was −0.22% a⁻¹ between the mid-1970s and ~2001 and was essentially unchanged at −0.23% a⁻¹ between ~2001 and ~2013. The most rapid shrinkage was measured in the southwestern sub-basin 5Z3 (Zhari Namco), at −0.33% a⁻¹ between the mid-1970s and ~2001 and −0.50% a⁻¹ between ~2001 and ~2013.

4. DISCUSSION

Our aim was to provide, for the first time over a large part of HMA, a multi-temporal inventory capable of shedding light not just on glacier shrinkage but on possible changes in rates of shrinkage. To achieve this aim in a reasonable time it was necessary to make some compromises. We found that automatic glacier mapping methods, such as band-ratio or normalized-difference methods, tend to misinterpret ice in deep mountain shadow as non-glacier (e.g. Fig. 3). Although the amount of misinterpretation varies with the choice of thresholds and other methodological details, we found that the advantage in productivity offered by automated methods is largely offset by the need for post-processing. Manual delineation of glacier outlines, time-consuming as it is, yields better results in complex terrain.

For practical reasons we decided not to attempt to map debris-covered ice, which implies that our estimates of total glacierized area are low. However, CGI-2 gives rather low estimates of debris cover (Table 6). Apart from basin 5Y6, where it is 7.8%, it exceeds 2.0% only in basins 5O and 5Q, and over the TP as a whole it is 3.4%. Thus our underestimates of total glacierized area are likely to be small in most basins. In some studies, (e.g. Shukla and others, 2009; Bajracharya and others, 2014), it has been reported that the proportion of debris-covered ice has been increasing at the expense of debris-free ice, and so our estimates of rates (as opposed to amounts) of total-area shrinkage may be exaggerated. However, Bajracharya and others (2014) showed for Bhutan that shrinkage rates for debris-free glaciers were indistinguishable from or moderately greater than those for all glaciers during three periods between 1977/78–1990 and 2010. Thus the effect of relative growth of debris cover on shrinkage rates at regional scales appears to be an open question.

We also decided not to subdivide glacier complexes into single glaciers in this study, which has no effect on the area change rates calculated at the basin scale. However, the three TPG datasets can be revised to develop a time series of glacier-by-glacier inventories since the 1970s. This work is now in progress for different mountain ranges.

The spatial pattern of shrinkage across the TP is illustrated in Figure 10. Negligible and moderate shrinkage rates calculated from TPG are located mainly in the northwest TP, where glaciers in 5Y6 (the southern Pamir and Hindu Kush eastwards to the Karakoram and Kunlun) showed the smallest change rate, −0.12% a⁻¹, and glaciers in 5Y5 in the northern TP changed by −0.17% a⁻¹. These rates are comparable with the insignificant rate reported by Holzer and others (2015) in part of the eastern Pamir. More rapid shrinkage rates were found in the southeast TP, for example glaciers in 5L presented the most rapid shrinkage rate, −1.60% a⁻¹ since 1994 and −0.43% a⁻¹ since the 1970s. The latter is statistically indistinguishable from the re-compiled result across HMA of Cogley (2016). In the IB from the 1970s to 2013, glaciers in 5Z1 and 5Z5 in the northern TP changed by −0.14% a⁻¹ and −0.17% a⁻¹ respectively, which is consistent with the results of −0.12 to −0.17% a⁻¹ from research by Wei and others (2014). Glaciers in 5Z4 and 5Z6 in the middle part of the TP experienced a moderate shrinkage by −0.20% a⁻¹. Glaciers in the southern TP, 5Z2 and 5Z3, reduced at the most rapid rates, by −0.24% a⁻¹ and −0.37% a⁻¹. Although glaciers in 5Z3 also showed rapid shrinkage rates in the present study, they were much slower than the result of −0.72% a⁻¹ obtained by Wei and others (2014).

Shrinkage rates clearly show a decrease from the southeast to the northwest of the TP. This spatial pattern is broadly consistent with patterns of area change reported earlier by Yao and others (2012), Li and others (2008) and Ding and others (2006) on the basis of smaller samples than that reported here. Interestingly, our broad spatial pattern of area change also agrees well with patterns of mass change reported by Yao and others (2012) on the basis of in-situ measurements and by Kääb and others (2015), Neckel and others (2014) and Gardner and others (2013) on the basis of altimeter measurements. The magnitudes of change, as opposed to the patterns, do not agree well in detail everywhere, but it is reasonable to attribute this to strongly variable glacial responses to variable climatic forcing at the mountain-range scale, as well as to the measurement uncertainties and methodological differences discussed next.

Obviously snow and clouds could have affected glacier delineation, especially in southeastern Tibet where cloudy days are rather frequent in every season. Indeed, cloudiness in southeastern Tibet is the main reason for our having to subdivide the M1970s epoch into three sub-epochs (see Tables 2–4).

GAMDAM has less glacierized area than TGP 2001 over the entire TP (Figs 11a, b). The time spans of the two inventories are similar, which means that differences in survey dates are not likely to explain the area differences. Nor is it likely that the area differences are due to different treatment of debris cover, which is omitted in TGP2001 and included in GAMDAM (Fig. 11b). For example, if we reduce the total glacierized area in GAMDAM by the 3.4% of TP-average debris cover from CGI2 (Table 6), the area difference (TGP2001 minus GAMDAM) becomes 14.4% instead of 10.4%, while in the basin with the most extensive debris cover, 5Y6, the difference increases from 3.6% to 8.5%. In fact, Figure 11 suggests clearly that many more ice-covered steep headwalls were excluded in GAMDAM than in TGP2001 (Fig. 11a). Nuimura and others (2015) acknowledged that ‘our exclusion of steep headwalls that are unaffected by glacier mass balance potentially discounts glaciers located on steep ground, resulting in an underestimation of total ice volume and median elevations in the GGI’. Arendt and others (2015) noted that the GAMDAM inventory excludes thin ice on headwalls and tends to have areas smaller than those measured in conformance with GLIMS guidelines. Moreover, fewer shaded glacier areas identified in GAMDAM (Nuimura and others, 2015) might also be one of the reasons for the glacier area difference. According to the GLIMS guidelines, perennial snow masses should be included as glaciers (Raup and others, 2007; Racoviteanu and others, 2009). The reasons why our inventory has 10% greater ice cover than GAMDAM require more thorough investigation, especially in the southeastern TP, but that is beyond the scope of this paper. Further glacier investigations relying on high-quality remote sensing data (optical or microwave) based on multiple platforms or in-situ investigations would be helpful.

Fig. 11. Glacier outline comparisons between the CGI2, RGI5.0 and TPG2001 datasets in the Hengduan Shan near 93°30′E, 30°20′N, where RGI 5.0 represents GAMDAM (GGI). (a) Ice outlines where TPG2001 is present while GGI is absent, mostly on steep headwalls. (b) Ice outlines where the GGI is present and TPG2001 is absent, mostly at the termini of debris-covered glaciers.

5. CONCLUSIONS

This study presented data on changing glacier extent on the Tibetan Plateau for the first time based on three epochs of ortho-calibrated Landsat satellite data from the 1970s to ~2013.

The comparative analysis of three mosaics shows that glacier area loss across the entire Tibetan Plateau is unequivocal. The glaciers in the M1970-76 subset showed shrinkage at −0.19% a⁻¹ during 1976–2001 and −0.20% a⁻¹ during 2001–13. Glaciers in the M1970-88 subset shrank at −0.34% a⁻¹ during 1988–2013, and in the M1970-94 subset at −0.42% a⁻¹ during 1994–2013. The most dramatic glacier shrinkage occurred in the Mekong basin (region 5L) in the southeast, with change rates for M1970s-76 glaciers of −0.43% a⁻¹ during 1976–2013 (−0.31% a⁻¹ during 1976–2001 and −0.73% a⁻¹ during 2001–13). The M1970s-88 glaciers in 5L changed by −0.78% a⁻¹ during 1988–2013, and M1970s-94 glaciers by −1.60% a⁻¹ during 1994–2013. Glaciers in the western and northern Tibetan Plateau changed at the slowest rates: −0.12% a⁻¹ in 5Y6, −0.14% a⁻¹ in 5Z1, −0.17% a⁻¹ in 5Y5 and −0.17% a⁻¹ in 5Z5 from 1976 to 2013. Overall, glaciers in externally drained catchments showed more dramatic shrinkage than those in the internally drained catchments. This shows a spatial gradient across the plateau of glacier shrinkage rate, decreasing from the southeast to the north and the west.

Comparisons of our time series with three earlier glacier inventories of High Mountain Asia indicated that basin-scale differences in estimates of glacier extent at comparable times were mostly moderate, exceeding ±10% in only one basin for TGP2001 minus CGI2 and five basins for TGP2001 minus GGI; for the whole Tibetan Plateau the respective differences were −1% and +10%. The differences might be due in part to the exclusion of debris-covered areas and perennial snow masses in our inventory, but the treatment of steep headwalls and inconsistencies between dates of imagery also play a role. The reasons for these discrepancies between large-scale inventories, discussed earlier by Paul and others (2013) among others for smaller glacier extents, merit detailed further investigation.