Inventory of Glaciers in the Shaksgam Valley of the Chinese Karakoram Mountains, 1970–2014

Abstract

The Shaksgam Valley, located on the north side of the Karakoram Mountains of western China, is situated in the transition zone between the Indian monsoon system and dry arid climate zones. Previous studies have reported abnormal behaviors of the glaciers in this region compared to the global trend of glacier retreat, so the region is of special interest for glacier-climatological studies. For this purpose, long-term monitoring of glaciers in this region is necessary to obtain a better understanding of the relationships between glacier changes and local climate variations. However, accurate historical and up-to-date glacier inventory data for the region are currently unavailable. For this reason, this study conducted glacier inventories for the years 1970, 1980, 1990, 2000 and 2014 (i.e., a ~10-year interval) using multi-temporal remote sensing imagery. The remote sensing data used included Corona KH-4A/B (1965–1971), Hexagon KH-9 (1980), Landsat Thematic Mapper (TM) (1990/1993), Landsat Enhanced Thematic Mapper Plus (ETM+) (2000/2001), and Landsat Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) (2014/2015) multispectral satellite images, as well as digital elevation models (DEMs) from the Shuttle Radar Topography Mission (SRTM), DEMs generated from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images (2005–2014), and Advanced Land Observing Satellite (ALOS) World 3D 30 m mesh (AW3D30). In the year 2014, a total of 173 glaciers (including 121 debris-free glaciers) (>0.5 km²), covering an area of 1478 ± 34 km² (area of debris-free glaciers: 295 ± 7 km²) were mapped. The multi-temporal glacier inventory results indicated that total glacier area change between 1970–2014 was not significant. However, individual glacier changes showed significant variability. Comparisons of the changes in glacier terminus position indicated that 55 (32 debris-covered) glaciers experienced significant advances (~40–1400 m) between 1970–2014, and 74 (32 debris-covered) glaciers experienced significant advances (~40–1400 m) during the most recent period (2000–2014). Notably, small glaciers showed higher sensitivity to climate changes, and the glaciers located in the western part of the study site were exhibiting glacier area expansion compared to other parts of the Shaksgam Valley. Finally, regression analyses indicated that topographic parameters were not the main driver of glacier changes. On the contrary, local climate variability could explain the complex behavior of glaciers in this region.

1. Introduction

Glaciers are key indicators of climate change and provide important freshwater resources [1]. Glacier changes (e.g., advance or retreat) have a significant impact on physical (e.g., local and/or global sea levels), biological (e.g., global ocean circulation and marine ecosystems) and social (e.g., water resources and tourism) systems [1]. The Shaksgam Valley is located on the northern slope of the Karakoram mountain region and situated in the transition zone between the Indian monsoon and the Central Asia dry arid climate zones. Investigations of glacier changes over times in this area are particularly useful to understand the local climate variation because of the sparse climate records (especially at high altitudes) [2]. Meltwater from the Shaksgam glaciers helps sustain the people and natural ecosystems in the region, contributes to hydroelectric power generation, and provides important water resources for the Yarkant river, one of the tributaries of the Tarim river (main artery for the oases at the northern margin of the Taklimakan desert) [3]. On the other hand, the Shaksgam glaciers are also one of the largest sources of glacial lake outburst floods in this region [4,5]. Interestingly, glaciers in this region have been reported to be quite stable (even advancing) compared to the global trend of glacier retreat, which also known as the “Karakoram Anomaly” [6]. Therefore, long-term and frequent monitoring of the Shaksgam glaciers is crucial to try to predict their future behaviors (and future water supply) [2,7]. However, field investigation of the spatial extent of these glaciers is often difficult to perform because the glaciers are located in extremely cold regions or at high elevations that are inaccessible or inhospitable to humans. Furthermore, some of the glacier extents are large and change slowly or advance quickly (some of them have surge activity) [6,8,9]. Hence, repeated measurements are needed across the entire region and over a long time period. Thus, multi-temporal and up to date glacier inventories are necessary for accurate estimation of glacier mass balance and also for hydrological modelling analysis [7,10,11,12]. Fortunately, satellite remote sensing data and image analysis techniques provide the means for delineating the areal extent of glaciers in time and space over large areas [12].

Until now, except for publicly available glacier inventories such as the First Chinese Glacier Inventory (FCGI) [13], the Second Chinese Glacier Inventory (SCGI) [14], Glaciers_cci [15] and the Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM) glacier inventory [16], only two studies have investigated the glacier changes in the Karakoram region, which encompasses the Shaksgam Valley [8,17]. However, there are still gaps from these existing studies. For example, the FCGI was done using aerial photography (acquired from the 1970s) that contained extensive fresh snow cover, which hampered the detection of the actual glacier boundaries and resulted in overestimation of some of the glacier extents [13,18] (see in Figure 1). In the Rankle et al. [17] study, in order to generated the updated glacier inventory (this inventory restricted to glaciers at least 3 km in length and more than 0.15 km² in area), glacier outlines from RGI (Randolph Glacier Inventory) were manually improved using Landsat images acquired from 2009 to 2011. In the case of Bhambri et al. [8] study, recent glacier inventory (only surge-type glacier was inventoried) was generated using same method with Rankle et al. [17] study based on the most recent Landsat images (2013–2015) and high resolution of Google Earth images. However, for the past period, both studies were only mapped the glacier terminus position changes (not the entire glacier outline) using multi-mission Landsat images and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images [8,17]. In addition, rather than performing the systematic image co-registration, both studies assumed the horizontal shift between a pair of satellite images (half pixels to 2 pixels). Importantly, only 17–54% of the glacier (compared to total number of glacier in this study) was inventoried in these studies (Figure 2). Thus, there is still a lack of glacier extent information in the historical time, and such inventory was unavailable for the Shaksgam glacier in the 1980/1990 periods and glacier outlines in 1970s from FCGI had larger uncertainty. Moreover, glacier extent information was still missing in the recent period, due to updated glacier inventory [8,17] not mapping the glacier in the entire Shaksgam valley. On the other hand, the Glaciers_cci, GAMDAM and SCGI inventories were performed using Landsat Thematic Mapper (TM)/Landsat Enhanced Thematic Mapper plus (ETM+) images acquired around 2000/2001 and 2009. Hence, variations in the quality (and biases) of glacier boundaries in these different data sets may make them unsuitable for performing a detailed analysis of glacier variations over time [19]. Consequently, a regularly updated and more systematic/consistent glacier area change assessment of the area is needed. To overcome the major shortcomings of the existing studies and data sets, this study performs a multi-decade glacier inventory for the 1970, 1980, 1990, 2000 and 2014 time periods using satellite remote sensing data (multispectral imagery and digital elevation data). Finally, glacier changes and inter-glacier relationships were estimated, analyzed and discussed for the 1970–2014 periods based on a set of topographic/climate factors.

2. Study Area and Datasets

2.1. Study Area

The Shaksgam Valley (Figure 3) is located at the northern slope of the Chinese Karakoram mountain range in the south-west of the Xinjiang Uyghur Autonomous Region, China. It has an area of 4644 km² and ranges in elevation from 3200 m to 8611 m (second highest peak in the world, K2). Due to the special location of the region, the Shaksgam Valley is under the complex influences of both the Indian monsoon and dry arid climate systems [20]. The westerlies (an important moisture source in winter) and advection of the moist air masses from the India monsoon are the main source of precipitation in this region [21]. The formation of large glaciers is a result of this special climate system (Figure 3). Therefore, the Shaksgam Valley has one of highest concentrations of large glaciers in mainland Asia [13]: e.g., eight glaciers are more than 50 km in length and more than 20 glaciers are at least 30 km long [13]. Debris-covered glaciers are also common in the Shaksgam valley, with most of the larger glaciers in this region being covered by thick debris [13]. However, there is no information about the exact debris thickness of these glaciers because of the inaccessible glacier terrain, which has resulted in almost no field observations [13,22]. The glaciers are an important contributor to the Indus and Yarkand rivers, and supports the livelihoods of around 130 million people [23]. Previous studies reported that Karakoram glaciers have stable or advancing terminus positions and surging behavior compared to the worldwide retreat of many mountain glaciers [8,17,24,25].

2.2. Datasets

Cloud-free satellite images with minimal seasonal snow cover, acquired at the end of the ablation season, were used for glacier mapping in this study (see

Table 1. Data used in this study.

Satellite	Scenes ID	Date (dd/mm/yy)	Resolution	Spectral Band	Utilization
Corona KH-4A	DS1025-1039DA015	08/10/1965	3.5 m	1 panchromatic band (PAN)	Glacier inventory for 1970
Corona KH-4A	DS1044-1023DA019	04/11/1967
Corona KH-4B	DS1107-1104DA005	30/07/1969
Corona KH-4B	DS1107-1104DA003	30/07/1969
Corona KH-4B	DS1107-1104DA002	30/07/1969
Corona KH-4B	DS1107-1104DA004	30/07/1969
Corona KH-4B	DS1107-1104DA003	30/07/1969
Corona KH-4B	DS1107-1104DA005	30/07/1969
Corona KH-4B	DS1115-1088DF201	16/09/1971
Corona KH-4B	DS1115-1088DF202	16/09/1971
Corona KH-4B	DS1115-1088DF203	16/09/1971
Corona KH-4B	DS1115-1088DF204	16/09/1971
Hexagon KH-9	DZB1216-500018L009001	22/07/1980	7.5 m	1 PAN	Glacier inventory for 1980
Hexagon KH-9	DZB1216-500361L006001	16/09/1980
Landsat TM	ETP148R35_5T19900629	29/06/1990	28.5 m 120 m	3 visible (VIS), 1 near-infrared (NIR), 2 shortwave infrared (SWIR), 1 thermal infrared (TIR)	Glacier inventory for 1990
Landsat TM	LT51480351993188ISP00	07/07/1993	30 m 120 m	3 VIS, 1 NIR, 2 SWIR, 1 TIR	Additional information for glacier identification (1990)
SRTM	SRTM3N35E075V2 SRTM3N35E076V2 SRTM3N36E075V2 SRTM3N36E076V2 SRTM3N36E077V2 SRTM3N35E077V2	02/2000	30 m		Additional information (Slope) for glacier identification (1990) Additional information (slope, plan curvature and profile curvature) for glacier identification (2000)
Landsat ETM+	LE71480352000168SGS01	16/06/2000	15 m 30 m 60 m	1 PAN 3 VIS, 1 NIR, 2 SWIR, 1 TIR	Glacier inventory for 2000
Landsat ETM+	LE71480352001202SGS00	21/07/2001	15 m 30 m 60 m	1 PAN 3 VIS, 1 NIR, 2 SWIR, 1 TIR	Additional information for glacier identification (2000)
Landsat OLI/TIRS	LC81480352014166LGN00	15/06/2014	15 m 30 m 100 m	1 PAN 3 VIS, 1 NIR, 2 SWIR, 1 TIR	Glacier inventory for 2014
Landsat OLI/TIRS	LC81480352015185LGN00	04/07/2015	15 m 30 m 100 m	1 PAN 3 VIS, 1 NIR, 2 SWIR, 1 TIR	Additional information for glacier identification (2014)
Terra ASTER	148-101-5_051020_L1A	20/10/2005	15 m 30 m 90 m	2 VIS, 2 NIR 6 SWIR 4 TIR	DEM generation and extracting additional information (slope, plan curvature and profile curvature) for glacier identification (2014)
Terra ASTER	148-101-4_110428_L1A	28/04/2011
Terra ASTER	149-100-7_110505_L1A	05/05/2011
Terra ASTER	147-102-2_110608_L1A	08/06/2011
Terra ASTER	148-101-5_120703_L1A	03/07/2012
Terra ASTER	148-100-2_140725_L1A	25/07/2014
Terra ASTER	148-101-2_140725_L1A	25/07/2014
AW3D30	N035E076 N035E077 N036E075 N036E076	2006–2011	30 m		Glacier parameters extraction
Soviet Military Topo maps	3 maps 7 maps	1975 1982	1:100,000		Comparisons for glacier inventory 1970–1980
GAMDAM Glacier Inventory		2000	30 m		Comparisons for glacier inventory 2000
Glaciers_cci		2000	30 m		Comparisons for glacier inventory 2000
Google Earth™ Images		2014	2–5 m		Comparisons for glacier inventory 2014

for the image data information). A detailed description of the datasets is provided in the following sections.

2.2.1. Corona KH-4A/B

For accurate delineation of the glaciers in the 1970s/1980s, declassified images from United States intelligence satellite missions such as Corona KH-4A/B and Hexagon KH-9 images were selected because of their high spatial resolutions (3–7 m) [18]. Corona KH-4AB contains a dual panoramic camera system (one camera looking forward and the other looking backward) with a separation angle of 30°, and offered the best ground resolution of about 1.8 m at nadir [26] (https://lta.cr.usgs.gov/declass_1). High performance photogrammetric film scanners were used (scanning process does not apply geo-corrections) to create digital products of Corona images in four segments at 7-micron (3600 dpi) with photogrammetric quality, and the image files are stored in TIFF format (https://lta.cr.usgs.gov/declass_1). For this study, two Corona KH-4A and 10 Corona KH-4B (acquisition time: 1965–1971, Table 1) images were used for generation of a glacier inventory for the year 1970.

2.2.2. Hexagon KH-9

The KH-9 mapping camera had a similar design to the National Aeronautics and Space Administration (NASA) Large Format Camera (LFC), with a 23 × 46 cm frame format camera and 30.5 cm focal length [27]. There are crosses (reseau points) on the Hexagon images, 23 in the x-direction and 47 in the y-direction. Hexagon KH-9 images provide a footprint of ~250 × 125 km and have approximately about 6–9 m spatial resolution [27,28]. Hexagon KH-9 photographs were scanned without any geo-corrections in two segments at 7 microns (3600 dpi) and available in an 8-bit TIFF file format (https://lta.cr.usgs.gov/declass_2). For this study, two Hexagon KH-9 images were utilized to produce glacier inventory around 1980.

2.2.3. Landsat Series Scenes, ASTER and DEMs

Orthorectified Landsat series scenes (TM, ETM+ and Operational Land Imager and Thermal Infrared Sensor (OLI/TIRS)) (Table 1) and digital elevation models (DEMs) from the Shuttle Radar Topography Mission (SRTM) (used for glacier mapping in 1990/2000), as well as mosaicked DEMs generated from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) stereo images acquired between 2005–2014 (representative the glacier geomorphometric parameters in 2014) were used for generation of glacier inventories in 1990, 2000 and 2014. Additional Landsat scenes from 1993 (Landsat TM), 2001 (Landsat ETM+) and 2015 (Landsat OLI/TIRS) (Table 1) were also used as complementary images during the digitization of the glacier outlines and to help differentiate seasonal fresh snow from permanent snow patches. 30 m spatial resolution SRTM and mosaicked ASTER DEMs were used to derive the geomorphometric parameters (GMP) such as slope, plan curvature and profile curvature, that are helpful for debris- covered glacier mapping [29,30]. Landsat TM, Landsat ETM+, and Landsat 8 images have six bands with 30 m spatial resolution in the visible (VIS), near-infrared (NIR) to shortwave infrared (SWIR) regions. Landsat ETM+ has a higher (60 m) spatial resolution thermal infrared (TIR) band compared to the Landsat TM TIR band (120 m) and Landsat OLI/TIRS TIR band (100 m). In addition, Landsat ETM+ and Landsat OLI/TIRS have a panchromatic band (PAN) with 15 m resolution that can be used to improve the glacier outlines derived from the semi-automatic processing technique.

2.2.4. DEM Data Used for Extracting Glacier Topographic Parameters

Glacier topographic parameters (e.g., minimum, maximum and median elevation, aspect, slope) have been found to be extremely useful for modelling, detection and change assessment of glaciers [31,32,33]. These topographic parameters have been extracted from the freely available SRTM or ASTER Global DEMs in past research. However, for this study we used the recently released (and also freely available) ALOS World 3D 30 m mesh (AW3D30) (http://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm) for the calculation of glacier parameters because it has been reported to have higher vertical accuracy (4.40 m root mean square error (RMSE) height accuracy) compared to the SRTM and ASTER Global DEM products [34,35]. The AW3D30 was completed by Japan Aerospace Exploration Agency (JAXA) collaborating with Nippon Telegraph and Telephone (NTT) DATA Corporation recently, and the 30 m product was derived by resampling the 5 m resolution (not freely available) global DEM dataset generated from the high resolution Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM), which was onboard the JAXA Advanced Land Observing Satellite (ALOS) [36].

2.2.5. Climate Data

Cloud cover, temperature and precipitation data were used to assess the possible climate influences on glacier changes in the Shaksgam Valley. The percentage of cloud cover was obtained from a Global 1-km Cloud Cover observation database which is derived from 15 years of twice-daily Moderate Resolution Imaging Spectroradiometer (MODIS) imagery at 1-km resolution [37]. Mean annual temperature and precipitation were derived from the China Meteorological Forcing Monthly Dataset (CFMD). CFMD was produced by merging the China Meteorological Administration stations data, Tropical Rainfall Measuring Mission (TRMM) satellite precipitation analysis data (3B42), Global Energy and Water Cycle Experiment-Surface Radiation Budget (GEWEX-SRB) downward shortwave radiation data, Princeton forcing data, and Global Land Data Assimilation System (GLDAS) data. This dataset currently covers the period 1979–2010 and has a 0.1-degree spatial resolution [38,39].

2.2.6. Validation Datasets

The evaluation of glacier outlines was performed based on Soviet Military Topo maps (from 1975 and 1982), Glaciers_cci, GAMDAM glacier inventory data, and freely available high resolution Google Earth™ images. The georeferenced Soviet Military Topo maps (1:100,000 scale), which show the glacier covered areas in 1975 and 1982, were used to validate the glacier outlines in 1970 and 1980. The GAMDAM glacier inventory data was used for comparison with the glacier outlines delineated for the year 2000 in this study. The GAMDAM glacier inventory was produced through manual delineation using 356 Landsat ETM+ scenes from the period 1999–2003 (images from 2000/2001 were used for the study area), in conjunction with a DEM and high-resolution Google Earth™ images [16]. Glaciers_cci was generated using conventional band ratio (band3/band5) with manual correction of debris cover, seasonal snow, shadow, water based on the Landsat TM/ETM+ images from the year 1998–2002 [15]. Furthermore, for evaluation glacier outline in 2014, Google Earth™ images from similar acquisition dates were used. Google Earth™ images provide multi-temporal higher spatial solution images and corresponding DEM.

2.2.7. Source of Datasets

Corona, Hexagon, Landsat series images and SRTM data were downloaded from the USGS Earth Explorer (EE) (USGS EE; http://earthexplorer.usgs.gov/). ASTER images (with minimal cloud cover and snow cover) were downloaded from the Japan Space Systems Earth Remote Sensing Division (ERSDAC) (http://gds.ersdac.jspacesystems.or.jp). AW3D30 was obtained from the JAXA Earth Observation Research Center (EORC) (http://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm). Soviet Military Topo maps were purchased from “Topographic Maps of the World” (https://mapstor.com/). Glaciers_cci was downloaded from the ESA CCI glaciers website (http://www.esa-glaciers-cci.org/). GAMDAM Glacier Inventory data was downloaded from Nagoya University (http://gamdam.cryoscience.net/). Global 1-km Cloud Cover observation database was available at http://www.earthenv.org/cloud. The CFMD dataset can be obtained freely from http://westdc.westgis.ac.cn/data/7a35329c-c53f-4267-aa07-e0037d913a21.

3. Method for Generation of Multi-Temporal Glacier Inventories

The generation of glacier inventories in this study follows the definitions developed within the GLIMS (Global Land Ice Measurements from Space) project. These definitions are adapted for the purpose of satellite-based glacier mapping based on official documents from the UNESCO (United Nations Educational, Scientific and Cultural Organization) guidelines for the compilation of the World Glacier Inventory [40]. For accurate estimation of glacier area change, only glaciers with areas larger than 0.5 km² were inventoried in this study because glacier outlines derived from Landsat images (the main source of glacier inventories in 1990, 2000 and 2014) may be more than 25% inaccurate for these very small glaciers due to the Landsat image pixel size [41]. The overall methodology for the generation of multi-temporal glacier inventories is shown in Figure 4, and the detailed processing steps are explained in the next Sections.

3.1. Glacier Inventory 1970/1980

3.1.1. Orthorectification of Corona Images

Information about the Frame Ephemeris Camera, Orbital Data camera, and Spacecraft parameters are classified for Corona missions [18]. Therefore, a traditional image georeferencing tool was used to orthorectify Corona images. At first, individual Corona image segments were merged and subsetted for the 15 subsets of Shaksgam Valley. Then, all subsets were co-registered by a spline projective transformation, which was performed based on ground-control points (GCPs) from a panchromatic band of Landsat ETM+ (2000) (performed in ESRI (Environmental Systems Research Institute, Inc., Redlands, CA, USA) ArcGIS 10.2). Similar features such as stable rock features, river junctions, and road stream junctions from Corona and Landsat ETM+ (2000) images were used for collection of GCPs, and more than 200 GCPs were collected for each subset image. The final orthorectified Corona images had a 3.5 m spatial resolution.

3.1.2. Orthorectification of Hexagon Images

KH-9 Hexagon images were treated as regular aerial photographs for orthorectification purposes. Thus, the Hexagon images could be orthorectified using the aerial imagery tool of COSI-Corr (http://www.tectonics.caltech.edu/slip_history/spot_coseis/index.html). Several steps were taken to orthorectify the Hexagon images. First, we made a mosaic of the two image halves of the Hexagon imagery. Second, we defined the interior orientation (IO) based on the reseau points on the Hexagon images. Third, we computed the exterior orientation (EO) based on the defined IO and selected GCPs (8 and 15 points) from a panchromatic band of Landsat ETM+ (2000) and SRTM. Fourth, orthorectified Hexagon images were produced by computing the mapping matrices using the optimized GCPs and EO. Finally, the orthorectified Hexagon images were co-registered again with the Landsat ETM+ (2000) imagery using a spline adjustment, and resampled to 7.5 m using ESRI ArcGIS. As a result, more than 250 GCPs were selected in the adjustment step.

3.1.3. Manual Delineation of Glaciers Based on the Corona and Hexagon Images

Manual digitization was the main method for generation of the glacier inventories in 1970 and 1980. Digitizing the debris-free glaciers was considerably easier because snow and ice present obvious brightness differences compared to the debris-covered glacier and surrounding objects. On the other hand, the boundaries of the debris-covered glaciers could be identified by their clearly recognizable geometric and topographic futures (ice crevasses, lateral and medial moraine) and their spectral differences (supraglacial ponds, exposed ice, river channel connect to the glacier terminus) from non-glacial areas in the Corona and Hexagon images (Figure 5).

3.2. Glacier Inventory 1990

For the circa 1990 period, glaciers were delineated from Landsat TM images and a binary slope map (<28°) derived from SRTM, following the method proposed by Alifu et al. [29]. This method maps the glaciers based on a combination of a band ratio image [TIR/(NIR/SWIR)] and slope characteristics. The main advantage of the [TIR/(NIR/SWIR)] band ratio method is that it can utilize the benefits of both optical and thermal infrared data sets. The (NIR/SWIR) part of the equation helps to identify clean glacier-ice, whereas the TIR band contributes to distinguishing supraglacial debris (debris cover on the glacier) from the surrounding periglacial debris (debris cover on the outside of the glacier margin) region based on the thermal differences. However, the use of [TIR/(NIR/SWIR)] band ratio alone can lead to several inaccuracies for the supraglacial debris classes in the areas where bedrock valley walls are located in the shade and/or in higher elevation areas. Therefore, additional information based on slope is included to eliminate these errors. Finally, manual editing was performed to improve the final glacier outlines.

3.3. Glacier Inventory 2000

Glaciers inventory for the year 2000 was mapped based on a combination of [TIR/(NIR/SWIR)] band ratio and result from the geomorphometric analysis using Landsat ETM+ (2000) and SRTM [30]. Unlike the method used for the glacier map in 1990, additional GMP such as slope (<28°), plan curvature, and profile curvature, were computed from the SRTM data and then, integrated using the iterative minimum distance statistical clustering technique implemented into the SAGA software (SAGAGIS: http://www.sagagis.uni-geottingen.de/html/index.php). Detailed steps of the methodology are presented in Alifu et al. [30]. Finally, manual improvements were conducted for final glacier outlines.

3.4. Glacier Inventory 2014

A similar approach to the glacier inventory generation (2000) was used for generating a glacier inventory in 2014 based on the Landsat 8 (2014) images and DEMs generated from ASTER (2005–2014) images. As explained previously, DEMs play a very important role in the accurate mapping of glaciers. However, the freely available global DEMs may not represent the current glacier terrain properties due to their less recent (or unknown) acquisition dates. For example, the SRTM DEM was generated from Space Shuttle Endeavour during the 11-day mission in February 2000. The other freely available global DEMs, ASTER Global DEM and AW3D30, were generated from stereo pairs images of unknown acquired dates, so, they also may not be representative of the glacier terrain properties in 2014 due to the changing topography of glaciers and the surrounding areas. Therefore, a more recent DEM was required for mapping glaciers circa 2014. Seven individual DEMs were generated from ASTER images using the SILCAST software, which can produce DEMs automatically without any ground control points (http://www.silc.co.jp/en/products.html) [42]. The individual raw DEMs were edited to remove the unnatural peaks and no data values (very few pixels). Then, corrected raw DEMs were mosaicked to produce a DEM of the whole study area using cubic resampling method. The final vertical accuracy of the mosaicked DEM was evaluated using comparisons of 109 points on the glacier-free terrain from ASTER GDEM, SRTM, and AW3D30. The evaluation results showed that DEM generated from this study has relatively good accuracy compared to other DEMs (

Table 2. Comparisons of 109 points from glacier-free terrain from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Digital Elevation Model (DEM) mosaic with reference DEMs.

DEMs	Total Points Number	Min	Max	Mean	DEMthis study RMSE
DEMthis study	109	3565	5360	4485	/
SRTM	109	3580	5356	4484	13
ASTER GDEM	109	3570	5377	4480	12
AW3DS30	109	3580	5336	4481	13

Note: Root mean square error: RMSE is estimated from the difference between DEM_{this study} and reference DEMs.

).

3.5. Complete Glacier Inventory

To complete the glacier inventory, the glacier-ice masses need to be separated into individual glaciers. Additionally, glacier parameters such as glacier area, debris-covered area, minimum (min) elevation, maximum (max) elevation, median elevation, mean slope angle, and mean aspect were calculated.

For this, first the glacier outlines were separated into individual glaciers based on manual inspection of the satellite images, the hillshade images derived from DEMs, and Google Earth™ high-resolution images. Through the ice-division processing, seasonal snow cover outside the glacier margins was also removed. Following the separation of individual glacier outlines, the debris-free glaciers area and supraglacial debris area of each individual glacier was calculated in ArcGIS (“Calculate Geometry” tool). The min elevation, max elevation, median elevation, and mean slope angle were also extracted on a glacier by glacier basis based on the AW3D30 data in ArcGIS (“Zonal Statistics” tool). The mean aspect or orientation of the glacier surface was calculated as the arctangent of the quotient of the sums of aspect sines and cosines of each of the glacier’s DEM grid cells [43,44]. Min elevation points were manually edited to the correct position and deleted if there were duplicate points (having the same elevation value). Then, in order to detect the glacier terminus change, spatial joins were performed in ArcGIS to calculate the distance between the two min elevations of points of each glacier in each time period.

4. Accuracy Assessment and Error Estimation

The potential errors of multi-temporal analysis come mainly from positional and mapping errors [45]. Co-registration of Landsat images was performed using ENVI software (Exelis Visual Information Solutions, USA) based on area-based matching method. The co-registration result showed a mean horizontal shift of about 0.35 pixel (10.5 m) for the Landsat TM images and, 0.25 pixel (7.5 m) for Landsat ETM+ and the Landsat OLI/TIR images. On the other hand, to assess positional accuracy for Corona and Hexagon images, 8–25 common location points such as river junctions and stable rocks were identified. Root mean square error (RMSE) of identified features (control points) was calculated by ArcMap (ESRI, USA) georeferencing tool using zero-order polynomial transformations (which is commonly used when the data is already georeferenced, but a small shift will better line up the data). RMSE of these features in each images showed that there was approximately 1.8 pixels (6.3 m) and 1.3 pixels (9.9 m) horizontal shift for Corona and Hexagon images. In this study, the uncertainty of glacier outlines was estimated by the buffer method recommended by Bolch et al. [45]. Hence, to estimate the uncertainty of glacier outlines, 6.3 m, 9.9 m, 10.5 m and 7.5 m buffer sizes were selected for the Corona, Hexagon, Landsat TM, and Landsat ETM and Landsat OLI/TIRS images, respectively. As a result, the mean uncertainty of the delineated glacier area was estimated as 2% for the Corona images, 3% for the Hexagon images, 4.6% for Landsat TM and 2.4%, 2.3% for the Landsat ETM+ and Landsat OLI/TIRS images (

Table 3. Error estimates for area (year), area change (period) and glacier terminus change (period).

Year	Errorarea Year (%)	Period	Errorarea Period (%)	Errorterminus Period (m)
1970	2.0	1970–1980	3.5	20.0
1980	2.9	1980–1990	4.3	42.0
1990	3.2	1990–2000	4.0	52.0
2000	2.4	2000–2014	3.3	50.0
2014	2.3	1970–2014	3.0	37.0

). The uncertainty of glacier area changes between each period was calculated from the uncertainty of the two differenced extents (Table 3) [46]. On the other hand, we estimated the errors in glacier terminus change (Table 3) based on the Equation (1), which proposed by Hall et al. [47]:

$$E=\sqrt{a^2+b^2}+e$$

(1)

where; E = error in glacier terminus change, a = pixel resolution of imagery a, b = pixel resolution of imagery b, e = horizontal shift (Note that significant glacier change is determined when rate of glacier changes is larger than estimated errors (Figures 8 and 12–14)).

In addition to positional errors, mapping errors of glacier outlines are a key source of uncertainty for glacier change monitoring. In our study, the degree of mapping error for each inventory was evaluated visually using the Soviet Military Topo maps (for the results from 1970 and 1980), the GAMDAM glacier inventory data (2000), and high resolution Google Earth™ images (~2014). Unfortunately, no appropriate reference data was found for visual validation of glacier outlines in 1990. Finally, visual comparisons show the glacier boundaries generated from this study were quite similar to reference datasets (Figure 6), and also our estimates of uncertainties were of a similar range to those of the previous studies [3,46,48].

5. Results

5.1. Updated Glacier (2014) Inventory and Glacier Characteristics

The glacier inventory in 2014 yielded 173 glaciers (including 121 of debris-free glaciers), covering an area of 1478 ± 34 km² (debris-free glacier covering an area of 295 ± 7 km²); 52 glaciers had debris cover on the ablation zones, covering an area of 128 ± 3 km². A similar number of glaciers were distributed in each sub-region (Figure 3b and Figure 7c). However, the eastern part (40% of total glacier area) of the Shaksgam Valley had the highest glacierized area, followed by the western part (37% of total glacier area) and the middle part (23% of total glacier area) (Figure 7c). About 8.6% of the total glacier area was covered by debris. Glaciers located in western part of the study area had the largest debris covered areas, followed by the middle and eastern parts. The elevations of glacier terminus positions ranged from 3972 to 5922 m, with an average of approximately 4960 m. The average median elevation of the glaciers was approximately 5578 m. The mean slope of the glaciers overall in 2014 was 26° (

Table 4. Topographic parameters of glaciers based on the glacier inventory 2014.

Inventory	Overall			Debris-Free Glacier			Debris-Covered Glacier
	Min	Max	Mean	Min	Max	Mean	Min	Max	Mean
Min elevation (m)	3972.0	5922.0	4960.0	4043.0	5922.0	5083.0	3972.0	5215.0	4676.0
Max elevation (m)	5560.0	7901.0	6168.0	5581.0	7207.0	6089.0	5560.0	7901.0	6351.0
Median elevation (m)	5013.0	6058.0	5579.0	5201.0	6058.0	5628.0	5013.0	5794.0	5464.0
Slope (°)	15.0	43.0	26.0	15.0	43.0	27.0	15.0	36.0	24.0
Area (km2)	0.5	344.6	8.5	0.5	84.7	2.4	0.6	344.6	22.8
Debris cover (km2)	-			-			0.03	41.5	2.5
Number of glacier	173			121			52
Total Area (km2)	1478			295			1184

).

In terms of the sizes of individual glaciers, 36% of glaciers were 0.5–1 km² in size, covering 3.1% of the total glacier area (Figure 7a). The highest number of glaciers (45%) were 1–5 km² in size, covering 11.5% of the total glacier area. Hence, 81% of glaciers were smaller than 5 km², but they covered less than 15% of the total glacier area (Figure 7a). On the other hand, the eight largest glaciers (greater than 50 km² and smaller than 400 km²) covered 65% of the total glacierized area. The distribution of glacier orientations shows that most of the glaciers faced north (N) (37%), north-east (NE) (21%), and north-west (NW) (15%), and they represented 43%, 39%, and 8% of total glacier area, respectively (Figure 7b).

On average, the size of debris-free glaciers was greater in the eastern part of the Shaksgam Valley than in the other sub-regions (Figure 7c). About 65% of the debris-free glaciers were 0.5–5 km² in size, and they were mostly orientated in the N and NE directions (Figure 7b). On the contrary, around 30% of the glaciers belonged to debris-covered glacier type, and they made up to 80% of total glacier area in the Shaksgam Valley (Figure 7a). Most of the debris-covered glaciers were oriented towards the NE, N and NW directions. In general, the glacier sizes as well as the range between the glacier min–max elevations were smaller for debris-free glaciers than for the debris-covered glaciers (Figure 6a and Table 4). On the other hand, the median elevation and terminus elevation of debris-free glaciers were higher than for debris-covered glaciers (Table 4). This demonstrated that the debris-covered glaciers typically had larger ablation areas than the debris-free glaciers. The larger glaciers (most of them belong to debris-covered glacier) in this region tended to expand to lower altitudes (longer length and relatively gentle slope), in contrast, smaller glaciers termini located at higher altitude (shorter length and relatively steeper slopes) (Figure 7a).

5.2. Glacier Area Change During Four Decades (1970–2014)

The results of the multi-temporal glacier inventory indicated that total area of glacier decreased by 22.25 ± 44.69 km² during 1970–2014, representing a loss of 1.48 ± 3.00% of glacier area (Figure 8 and Figure 9). In detail, results for the different sub-regions showed that glacier area in the middle and eastern parts of study area decreased by 0.82 ± 3.00% and 3.26 ± 3.00%, while the glacier covers in the western part of Shaksgam Valley exhibited very little change (+0.08 ± 3.00%) (Figure 8). However, the estimated area change was less than or similar to the estimated error (Table 3), and thus the glacier area changes of the overall region as well as sub-regions was not significant. In detail, Figure 8 illustrates that except for several glaciers located in the eastern part of study area, which showed significant area decreases compared to the western part of the study site (where the glacier areas mostly increased), most of the glaciers area were stable over this period.

Further analysis of individual glacier area changes by sub-region showed that debris-covered glaciers (>50 km²) were stable in size in all sub-regions. About 60% of shrinking glacier (area decrease: −5–~26%) in eastern part of study are belonged to small debris-covered glaciers and small debris-free glaciers (<2 km²). On the other hand, small glaciers (most of them <1 km² and ~26% of total glacier in middle part, Figure 8) in the south region of middle part of the study area, and medium-sized glaciers (bigger than 2 km² and less than 11 km², and ~27% of total glacier in the western part) in the western parts of Shaksgam Valley significantly increased in area, by 5–28%.

5.3. Decadal Changes in Glacier Area between 1970 and 2014

The overall results of the glacier inventories for the five different periods are shown in

Table 5. Comparison of glacier inventories from this study (glacier area: km²).

Comparison		TGN	Total GA	Mean GA	SDGA
1970	Overall	176	1501	8.5	30.8
	Debris-free glacier	123	329	2.7	8.2
	Debris-covered glacier	53	1171	22.1	52.3
1980	Overall	177	1507	8.5	30.9
	Debris-free glacier	123	327	2.7	8.2
	Debris-covered glacier	54	1180	21.9	52.2
1990	Overall	173	1451	8.4	30.5
	Debris-free glacier	120	291	2.4	7.9
	Debris-covered glacier	53	1159	21.9	51.4
2000	Overall	173	1446	8.4	30.8
	Debris-free glacier	120	289	2.4	7.8
	Debris-covered glacier	53	1157	21.8	51.9
2014	Overall	173	1478	8.5	31.4
	Debris-free glacier	121	295	2.4	7.7
	Debris-covered glacier	52	1184	22.8	53.5

Note: TGN, total glacier number; GA, glacier area; SDGA, standard deviation of glacier area.

, Figure 10 and Figure 11.

The largest glacier area changes (−3.75 ± 4.3%) occurred during the 1980–1990 period, while minor glacier area changes occurred over the 1970–1980 (+0.44 ± 3.5%) and 1990–2000 (−0.36 ± 4.0%) periods. In the most recent period (2000–2014), the total glacier area increased by 2.3 ± 3.0%. Decadal glacier area changes overall and in the different sub-regions indicated that no significant glacier area changes occurred. However, further analysis of individual glaciers revealed that most were stable in size during the 1970–1980 period except for small glacier (<2 km²) decreased in area which contain 35% of glacier in eastern part of study area (Figure 10a) and 1990–2000 periods (Figure 10c). The 1980–1990 period had the highest number of glaciers with significant decreases in area, and these shrinking glaciers were distributed in all sub-regions (Figure 10b). In the 2000–2014 period, the properties of individual glacier area changes showed that most glaciers were stable in size or even increasing in glacier area, with the exception being some very small glaciers that decreased in area (mostly located in the eastern sub-region) (Figure 10d).

More detailed analysis of individual glacier changes in the different sub-regions between 1970 and 2014 indicated that large numbers of glaciers decreased in area (significant area change) in the eastern part of the study area in almost all periods, whereas glaciers located in the western part of the study area were stable in size (some glaciers even showing increases in area) during 1970 to 2014 (Figure 10). On the other hand, analyzing glacier area changes for the different size classes of glaciers showed that small debris-covered and small debris-free glaciers area decreased in size between 1970 and 2000 (Figure 10). In the case of larger glaciers (>50 km²), they tended to be stable or having slight increases in glacier area from 1970 to 2014 (Figure 10).

5.4. Glacier Terminus Change

The distance between the two points representing the minimum elevation of glaciers in different time periods indicates the glacier terminus change. In the entire study area, 32% of glaciers advanced, 22% of glaciers were stable, and 46% of glacier experienced retreat during the 1970–2014 period (Figure 12). However, compared to the advancing glaciers, retreating glaciers tended to be smaller glaciers. Comparing the debris-free and debris-covered glaciers, the number of advancing and retreating debris-free glaciers were more than debris-covered glaciers. Nevertheless, debris-covered glaciers showed a larger total terminus advance (range from 70 to 1400 m) and a smaller retreat compared to debris-free glaciers (Figure 12). The terminus changes of glaciers in each sub-region indicated that the middle part was close to equilibrium between glacier advance and retreat (Figure 12b). Retreating glaciers dominated in the eastern part of Shaksgam Valley, and most of the advancing glaciers were located in the western part of the study area (Figure 12c).

For the 1970–1980 and 1990–2000 periods, 47% and 56% of glacier terminus were stable, respectively, and 21% and 18% of glaciers’ terminus advanced (Figure 13a–c and Figure 14a–c). Glacier terminus changes from 1970 to 2014 indicated that 42% of glaciers (more than half being small debris-free glaciers) retreated during the 1980–1990 (Figure 13d–f) period. Compared to other periods, glacier advance dominated in the recent period (2000–2014) (Figure 14d–f). The number and magnitude of retreating glaciers was greater (in number) and larger (in length) than advancing glaciers during the 1980–1990 period. In detail, retreating of debris-free glaciers occurred in almost all sub-regions and periods (except in the western part during 1980–1990, and the middle and western part during 2000–2014). Debris-free glaciers in the eastern part of the study area showed the largest retreat in all periods. In the case of debris-covered glaciers, the number of retreating debris-covered glaciers was less than for debris-free glaciers, while the total magnitude of advancing debris-covered glaciers was more than for debris-free glaciers. Debris-covered glaciers experienced the most significant advances in the middle (3 glaciers advanced more than 1 km) and western (3 glaciers advanced more than 1 km) parts of the Shaksgam Valley during the most recent decade (2000–2014). Most of the advancing glaciers were oriented in the N, NE, and NW directions with a mean slopes of 26° (90% of the glaciers were in the 14–31° slope range).

6. Discussion

6.1. Glacier Change

In general, overall glacier changes over time in the Shaksgam Valley were quite stable (especially during the 1970–1980 and 1990–2000 periods). Similar to most of the Himalayan glaciers [49], our observation shows that glaciers in Shaksgam Valley retreated in all sub-regions in the 1980–1990 period [50]. Most of the retreating glaciers in this period belonged to the small glacier class (debris-covered and debris-free), and some of them were tributaries of large debris-covered glaciers (i.e., Figure 15b, 1) where the main reason for large debris-covered glacier area decreased which is shown in Figure 10b (in green). The smaller glaciers are, therefore, more sensitive to climate changes due to their location as lower median elevations and smaller ablation regions. Nonetheless, several individual glaciers were still advancing during this period (i.e., Figure 15b, 3). Debris-free glaciers experienced much greater retreats than debris-covered glaciers (except for in the western part of the study area), possibly because the thick debris layer on the debris-covered glaciers, typically concentrated at lower elevations where most of the melting is expected to occur, slowed down the retreats of the debris-covered glacier. Interestingly, the glacier area increased in the most recent period. Several glaciers which belonged to the middle size classes (5–10 km²), mostly located in the western and middle parts of the Shaksgam valley, advanced and merged with large glaciers (i.e., Figure 15b, 2 and Figure 15c, 2). This trend may likely be driven by changes in topography, debris variations (unfortunately no information about debris thickness is available), and the strong effect of regional climatology.

The trend of glacier area changes in this study shows a relatively similar trend to the study by Rankl et al. [17] and other glacier elevation change studies [51,52]. Glacier elevation change studies from Rankl and Braun [51] (Figure 2 in [51]) and Zhou et al. [52] (Figure 4 in [52]) also indicated that glaciers in the western part of Shaksgam Valley have been nearly stable (even though some glaciers were thickened) compared to the other parts of study area. Comparisons of glacier outlines and glacier area (Rankl et al. [17] data was downloaded from https://doi.pangaea.de/10.1594/PANGAEA.873278 and Bhambri et al. [8] data was provided by the first author through a personal communication) from this study and the previous studies were shown in Figure 16. Comparisons of glacier outlines indicated that even using relatively higher resolution of ASTER/Landsat images (15–30 m), precisely detection of debris-covered glaciers boundary by manual digitization (strongly reliant upon spectral information) still faces difficulty when identifying the completely debris-covered glacier termini (Figure 16a). In addition to differences mentioned above (i.e., Figure 16a–c), Bhambri et al. study has overestimated glacier area in this region because the glacier outlines contain the large area of seasonal snow cover (Figure 16d–e). On the other hand, the Rankl et al. study was underestimated the glacier area due to some tributary glaciers were not digitized (Figure 16d). These are the main area differences from previous research compared to this study.

A recently published review study reported that only four area change studies are available from the Karakoram Range covering different regions [9]. However, only one work studied Shaksgam glaciers change, but that study [53] must be regarded with caution because estimation of the glacier area was based on topographic maps (which is source of FCGI)/Landsat images and its uncertainty was not assessed [9]. Therefore, the result of this study provided a comprehensive analysis of glacier area change and detailed multi-temporal glacier inventories compared to the previous study. Moreover, results from our study should have more accurate glacier change estimates than previous work due to (1) the use of high resolution declassified images for the 1970/1980 periods and, (2) the use of semi-automatic glacier mapping methods based on the spectral, thermal and landscape properties of glaciers for more recent periods. However, our study excluded very small glaciers (<0.5 km²) because those glaciers present larger uncertainty when delineated using Landsat images [40]. Finally, the multi-temporal glacier outlines (including topographic parameters) of the Shaksgam Valley will be published with this paper and will be made freely available.

6.2. Glacier Area Change with Glacier Topographic Parameters

Glacier changes can be strongly influenced by topographic characteristic and climate variations. Thus, multi-regression analysis was applied to examine the relationship between the magnitude of glacier change (percentage of glacier area change bigger than 5%) with various topographic parameters.

Regression analysis results are illustrated in

Table 6. Regression between glacier area change and glacier parameters (84 observations) and regression between glacier area change and debris-cover (52 observations).

Regression	R Square	Coefficient	p Value
Glacier area	0.087	0.042	0.590
Min elevation		−0.001	0.817
Max elevation		0.001	0.898
Median elevation		0.002	0.684
Slope		−0.134	0.196
Aspect		−0.003	0.929
Debris-cover	0.016	0.101	0.110

and

Table 7. Regression analysis between glacier terminus advance (56 observations, left) and retreat (77 observations, right) with glacier parameters.

Regression	R Square	Coefficient	p Value	R Square	Coefficient	p Value
Glacier area	0.172	−2.417	0.050	0.264	−2.103	0.414
Min elevation		−0.401	0.155		−0.725	2.542
Max elevation		−0.206	0.304		−0.177	0.211
Median elevation		0.296	0.462		0.465	0.058
Slope		−2.373	0.796		−6.110	0.252
Aspect		−0.202	0.586		−0.212	0.303

. According to the results, the p value for all variables were found to be not significant at <0.001, while the R² value of the regression value was also relatively low: R² = 0.087 (R² = 0.016 for debris cover). On the other hand, the result of regression analysis of glacier advance/retreat also showed a low R² and no significant p values. Overall, the regression results indicate that the topographic parameters were not the main driver of glacier change in the Shaksgam Valley.

6.3. Glacier Area Change with Climate Considerations

Climatic factors such as cloud cover, precipitation, and temperature may have also caused glacier changes in Shaksgam Valley. In general, incoming short-wave radiation is the major source of energy for glacier change [54]. However, cloud cover may alter the amount of incoming short-wave radiation as well as air temperatures [54]. The percentage of cloud cover in the Shaksgam Valley was obtained from the online Global 1-km Cloud Cover observation database [37]. From this dataset, we found that the mean annual percent cloud cover and mean monthly (May–October) percent cloud cover were higher in the western part of study area compared to the other regions (Figure 17). Furthermore, most of the areas with frequent cloud cover were located on the accumulation zones of the glacier. Previous studies based on the weather observation stations and global meteorological reanalysis data, also reported an increase of cloudiness in the Karakoram region [54,55]. Moreover, the previous research also reported that increases in cloudiness led to increases in precipitation, decreases in incoming solar radiation, and a reduction in the air temperatures [54,55].

Moreover, variation of cloudiness is concurrent with temperature and precipitation change data derived from the CFMD dataset during 1979–2015 (Figure 17 and Figure 18). Further analysis of the spatial distribution of mean annual temperature and precipitation from CFMD (Figure 18) indicated that lower temperature and higher precipitation were dominated in the western part and south of the middle part of the Shaksgam Valley compared to the eastern part. Spatial variations of mean annual precipitation also revealed increased precipitation in the eastern part of the study area in the recent period (Figure 18h). Similar results were reported based on observation stations [56] (such as Tashkorgan, which is the nearest meteorological station to the Shaksgam Valley, see Figure 1) and satellite-based estimations (TRMM) [57]. Based on the meteorological stations, Yaning et al. concluded that annual mean temperature slightly increased during 1960–2010 [56]. Moreover, the decadal mean and decadal summer mean temperature even showed decreases in the recent decade. On the other hand, the result of precipitation change derived from observations as well as satellite data based showed a steady increase from 1960–2010 [56,57,58]. Precipitation estimates from satellites also indicated that the amount of precipitation in the glacial region had been underestimated [57]. Bashir et al. further confirmed that similar in precipitation and temperature were observed in the upper Indus basin located closely to the study area [55].

Results revealed that most of the shrinking glaciers are small glaciers, and small glaciers are more sensitive to the climatic variations than larger glaciers [59,60,61], therefore, climate variations can explain the glacier area decreases dominant in the eastern part of the Shaksgam valley during the 1970 to 2000 period (Figure 8, Figure 10 and Figure 18). On the other hand, long-term higher precipitation and lower temperature interpreted the climate variations played an important role in glacier advances in the western part of study area [62]. Thus, the glacier changes in Shaksgam Valley probably have a strong relationship with specific regional climate variations (combination of temperature and precipitation).

6.4. Limitations of the Study

There are several limitations in the proposed method used in this study, which can be summarized as:

• Experience from transferability of this method to other regions [29,63] indicated that the [TIR/(NIR/SWIR)] band ratio image misclassified vegetation cover near the glacier tongue as supraglacial debris, particularly in the case of debris-covered glacier covered by vegetation.
• Rock glaciers and ice-cored moraines are not possible to detect by the proposed method due to their much smaller size and overall area being entirely covered by rocky materials compared to the debris-covered glacier. Importantly, the thermal properties of rock glaciers and ice-cored moraines are undetectable because of a very thick deposit layer [64].
• Small hanging debris-covered glaciers had the highest uncertainty in glacier delineation mainly due to the limited image resolution. Therefore, manual editing steps are still required for small debris-covered glaciers and river channels connected with a glacier tongue and isolated areas that are outside the main body of the glacier.

7. Conclusions

This research takes advantage of various freely-available remote sensing data sets to generate multi-temporal glacier inventories (1970, 1980, 1990, 2000 and 2014) of the Shaksgam Valley located in the northern slope of the Chinese Karakoram region. As of 2014, 173 glaciers (>0.5 km²) (including 121 debris-free glaciers), covering an area of 1478 ± 34 km² (debris-free glacier area: 295 ± 7 km²) were inventoried. About 10% of the total glacier area was debris-covered. Overall, the total glacier area change was stable from 1970–2014. Nonetheless, individual glacier changes showed significant variability. Glaciers in the Shaksgam Valley still showed significant advances of terminus positions, especially in the recent decade, compared to the worldwide trend of glacier retreat. Notably, glaciers located in the western part of the study area experienced significant levels of increase in glacier area. The results from a regression analysis indicated that glacier changes in the Shaksgam Valley were not strongly related to the local topography, and were more likely due to local climate variations (e.g., changes in cloud cover, precipitation, and air temperatures). Unfortunately, it was not possible to investigate how the debris thickness of glaciers affected the rates of glacier change in this region due to the lack of debris thickness data. However, we did find clear evidence that the variation of small glacier-ice extent was more sensitive to the climate variations compared to the larger glaciers. Finally, because (multi-temporal) glacier inventory data is a standard requirement for various climatological and hydrological studies, the glacier inventory data from this research has been made freely available as supplemental material to this paper.