A Surging Glacier Recognized by Remote Sensing on the Zangser Kangri Ice Field, Central Tibetan Plateau

Abstract

A glacier surge, which is quasi-periodic and involves rapid flow, is an abnormal glacier motion. Although some glaciers have been found to be surging, little is known about surging glaciers on the Tibetan Plateau (TP), especially the Central and Northern TP. Here, we found a surging glacier (GLIMS ID: G085885E34389N) on the Zangser Kangri ice field (ZK), Central TP, by means of the digital elevation models (DEMs) from the Shuttle Radar Topography Mission (SRTM), TanDEM-X 90 m, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs, and High Mountain Asia 8-m DEM (HMA), combined with Landsat images and the Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE) dataset. This surge event was confirmed by the crevasses, shear margin, and visible advancing snout shown in the Landsat images produced since 2014 and the HMA. The inter-comparison of these DEMs and the surface velocity changes showed that the surge event started between October 2012 and January 2014. The glacier may have also surged in the 1970s, based on a comparison between the topographical map and Landsat images. The glacier mass balance here has been slightly positive from 1999 onward (+0.03 ± 0.06 m w.e.a⁻¹ from 1999 to 2015, +0.02 ± 0.07 m w.e.a⁻¹ from 1999 to December 2011), which may indicate that the ZK is located on the southern edge of the mass balance anomaly on the TP. Combining with other surging glaciers on the Central and Northern TP, the relatively balanced mass condition, large size, and shallow slope can be associated with glacier surges on the Central and Northern TP.

1. Introduction

Surging glaciers show quasi-periodic rapid motion [1]. The cycle of a surge usually involves a short active phase with rapid movement and a long quiescent phase with slow velocity [2]. Glacier velocity is 10–1000 times higher in the active phase than in the quiescent phase. This rapid motion transfers a large volume of ice and snow downward, leading to thinning on the upper reservoir area and thickening on the lower receiving area, which is a characteristic of surges [3,4]. Moreover, there are other indicative characteristics, such as crevasses, distorted medial moraine, and terminus advance [5,6]. Two major types of surges, the Alaskan-type and the Svalbard-type, are thought to have distinctions. Under the hydrological mechanism, the transition of the drain system from an open basal-tunnel system to a closed basal-cavity system initiates an Alaskan-type surge [7]. Under the thermal mechanism, the warming of basal ice caused by geothermal flux, surface temperature, and thick ice can initiate a Svalbard-type surge [8,9]. The Alaskan-type surge generally starts in the winter and terminates in the summer with a short active phase (several months to 1–3 years) [10]. The Svalbard-type surge has a long active phase (>3 years) and its initiation shows no clear seasonality [11]. The velocity changes also differ: the deceleration period is shorter compared with the acceleration period for the Alaskan-type surge; this is the inverse for the Svalbard-type surge [11]. However, according to studies in the Karakoram and Alaska Range [12,13,14], glacier surges can show a spectrum of speedup and flow instability, demonstrating inexplicit distinctions and various controlled processes. The mass balance [15], glacier geometry [16,17], and the deformable bed [4,18] can also control the surge. Consequently, the surge mechanism needs to be better understood.

As of 2013, there were 2317 documented surging glaciers, which account for ~1% of global glaciers [19,20]. Surging glaciers mainly concentrate in Alaska/Western Canada, Arctic Canada [21], Svalbard [22], and Pamir–Karakoram in Central Asia [23,24]. Although representing a limited portion of global glaciers, surging glaciers play an important role in improving our understanding of glacier flow and evolution [25]. Surging glaciers may also produce some catastrophic effects, as the fast-moving ice can change the landform and even threaten inhabitants and herds at the downstream part [26,27]. The Tibetan Plateau (TP) is the most glacierized land area outside the polar regions [28]. Several surging glaciers were recorded on the Southeastern TP in the twentieth century [29]. In the vast Central and Northern TP, with the coldest and driest climate in the TP [30], glaciers are thought to be stable and are classified as extreme continent/subpolar glaciers with a frozen bed and limited movements [31]. The harsh environment and huge logistic costs lead to few in situ observations here [28]. However, remote sensing data can shed light on the glaciology study in this region. Some glaciers have shown surging characteristics in the Muztag, Western Kunlun Mountains (WKM), Puruogangri, Geladaindong, and Xinqingfeng and Malan ice caps (Figure 1a) since 2000 [17,32,33,34,35,36].

Zangser Kangri (ZK) is located on the Central TP (Figure 1a). From May to September, the temperature is positive and ZK receives more than 80% of its annual precipitation (Figure A1). A total of 54 glaciers circle around the main peak, their equilibrium-line altitudes (ELA) ranging from 5700 to 5940 m above sea level (ASL) according to the first and second Chinese glacier inventories [37,38]. They have suffered a slight area reduction rate of −0.15% a⁻¹ to −0.10% a⁻¹ since the 1970s compared with other glaciers on the TP [39,40]. The glacier mass balance here was +0.37 ± 0.25 m w.e.a⁻¹ during 2003–2009 [41] or nearly balanced since 2000 [42,43]. Glacier G085885E34389N (the Global Land Ice Measurements from Space (GLIMS) ID) is an eastward-flowing glacier with a length of 13.5 km and an area of 35.63 km² in 2006 in the second Chinese glacier inventory.

In this study, a surging event on glacier G085885E34389N was identified using sequential digital elevation models (DEMs) and Landsat images. Its characteristics were recognized and studied from multiple remote sensing data. The velocity changes of this surge event were derived using the Co-registration of Optically Sensed Images and Correlation (COSI-Corr) software package and collected from the Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE) dataset. Characteristics of glacier mass balance, size, and slope are associated with surges on the Central and Northern TP.

2. Materials and Methods

2.1. Data

2.1.1. SRTM DEM

We collected 1-arcsecond Shuttle Radar Topography Mission C-band DEM (SRTM C DEM) from the United States Geological Survey (USGS) and ~25 m Shuttle Radar Topography Mission X-band DEM (SRTM X DEM) from the Deutsches Zentrum für Luft- und Raumfahrt (DLR). These DEMs were acquired from 11 to 20 February 2000 simultaneously, but their swath widths are different: 225 km for C-band data and 50 km for X-band data. Thus, the SRTM X DEM provided less coverage compared to SRTM C DEM in our study area (Figure 1b). Elevations in SRTM C DEM were orthometric heights relative to the EGM96 geoid. In SRTM X DEM, elevations were ellipsoidal heights relative to the WGS84 ellipsoid, which meant we should convert one from another when comparing these two DEMs. These two DEMs were used for elevation change and glacier mass balance estimation.

2.1.2. TanDEM-X 90 m

TanDEM-X 90 m (TDM90) was downloaded from the Earth Observation Center (EOC) of DLR. This data is a free 90 m global DEM generated from the bistatic X-band interferometric synthetic aperture radar (SAR) system of the TerraSAR-X satellite (TSX) and the TanDEM-X satellite (TDX), launched in June 2007 and June 2010, respectively. Elevations in this data are integrated values of multiple TanDEM-X DEMs from December 2010 to January 2016, and they are ellipsoidal heights relative to the WGS84 ellipsoid. We checked the xml file of the TDM90 data used in this study and found that glaciers on ZK were covered by 9 TerraSAR-X/TanDEM-X raw tiles. They were acquired from August 2011 to October 2012 (

Table A1. Information on raw tiles in TDM90 covering the ZK glaciers.

ID	Acquisition Start Time
1025046_11_PO_ITP_14	2011-08-19T12:11:54.656099Z
1025046_12_PO_ITP_14	2011-08-19T12:12:01.656028Z
1042359_30_PO_ITP_06	2011-11-04T12:11:56.053254Z
1057021_11_PO_ITP_12	2012-03-04T12:11:52.092323Z
1059493_03_PO_ITP_11	2012-01-31T12:11:52.456143Z
1059493_04_PO_ITP_11	2012-01-31T12:11:59.182974Z
1081220_08_OP_ITP_16	2012-10-10T12:11:59.928249Z
1081220_09_OP_ITP_16	2012-10-10T12:12:06.928258Z
1083733_07_OP_ITP_16	2012-09-07T12:12:00.078820Z

). This DEM was used for elevation change and glacier mass balance estimation. TDM90 leads to a time uncertainty (~1 year) when estimating the mass balance.

2.1.3. High Mountain Asia 8-Meter DEM

High Mountain Asia 8-meter DEM (HMA) was downloaded from the National Snow and Ice Data Center (NSIDC). This DEM dataset was generated from level-1B (L1B) panchromatic stereo images over the High Mountain Asia of WorldView-1, WorldView-2, WorldView-3, GeoEye-1, and Quickbird-2. It was processed by the Ames Stereo Pipeline (ASP) from the National Aeronautics and Space Administration (NASA). The major temporal coverage of this dataset was from 2013 to 2016 [43,44]. Elevations in this dataset are relative to the WGS84 ellipsoid. The tile used in this study was acquired in October 2015. These data were used for elevation change and glacier mass balance estimation. There are some voids on glaciers in the selected HMA (Figure 1b).

2.1.4. ASTER DEMs

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images contain 14 bands with resolutions ranging from 15 to 90 m. Band 3 has two views, the nadir and the backward, which can be used to generate a DEM. ASTER DEMs (the AST14DMO product) were generated by the new Land Processes Distributed Active Archive Center (LP DAAC), NASA. However, some researchers showed that this product is affected by jitter with erroneous elevation values [36,45]. Therefore, two well-processed ASTER DEMs from [42] were used. They represent elevations in May 2014 and March 2015. These two ASTER DEMs were matched to SRTM C DEM, relative to the EGM96 geoid (personal communication with Fanny Brun). Given the large voids in these two ASTER DEMs (their coverages can be seen in Figure 2c,d), they were used for determining the start time of the surge.

2.1.5. Landsat Images and Topographic Map

We downloaded 1 Landsat multi-spectral scanner (MSS) image from 1977, 4 Landsat thematic mapper (TM) images from 1991 to 2006, 6 Landsat enhanced thematic mapper plus (ETM+) images from 1999 to 2012, and 10 Landsat operational land imager (OLI) images from 2013 to 2020 from USGS. All Landsat images were processed by radiometric and geometric corrections under the same Universal Transverse Mercator (UTM) 45N projection. Images from July to November and with less than 10% cloud cover were preferentially selected but we were obliged to work with images from other months that had more cloud coverage. These images were used for digitizing glacier outlines, delineating the flow line of glacier G085885E34389N, and measuring velocities. One scanned topographic map (scanned resolution of ~9 m [39], a scale of 1:100,000) was collected for flow line delineation of glacier G085885E34389N. It was generated from aerial stereo pairs acquired in December 1971 by the Survey and Mapping Bureau, General Staff Department of the People’s Liberation Army.

2.1.6. Global Land Ice Velocity Extraction from Landsat 8 (GoLIVE)

The GoLIVE dataset (Version 1.1) is a land ice velocity compilation that has been generated from pairs of Landsat 8 panchromatic images since May 2013, and it can be downloaded from the National Snow and Ice Data Center (NSIDC) [46]. Velocities were derived using the image correlation software, Python Correlation (PyCorr), from images with the same orbit, land ice >5 km² (according to the Randolph Glacier Inventory), and cloud cover lower than 50% [47]. This dataset provides velocities of most glaciers within the latitude range 82° S to 82° N from 2013 in Network Common Data Format (NetCDF) with a sample spacing of 300 m. Velocity data in the Landsat World Reference System-2 of 142/036 were used. These data are under the UTM 45N projection and have time intervals from 16 to 96 days. The variable vv_masked in the NetCDF, which represents the velocity with high confidence and high accuracy, was chosen. The details of all remote sensing data and their uses are listed in

Table A2. Information on remote sensing data.

Data	Product ID/Date	Pixel Size/Scale	Cloud Cover (%)	Purpose
SRTM C DEM	SRTM1N34E085V3	1-arcsecond		Elevation change
SRTM X DEM	E080N30	~25 m		Elevation change
HMA	HMA_DEM8m_AT_20151002_0649_1020010044807000_1020010043B63700	8 m		Elevation change
ASTER DEM	AST_L1A#00305052014050415_05052014194701	30 m		The start time of the surge
ASTER DEM	AST_L1A#00303122015051010_03132015141806	30 m		The start time of the surge
ASTER DEM	AST_L1A#00303122015051019_03132015141819	30 m		The start time of the surge
Topographical maps	9-45-52/December 1971	1:100,000		Flow line delineation
Landsat MSS	LM02_L1TP_152036_19770303_20180422_01_T2	60 m	1.00	Flow line delineation
Landsat TM	LT05_L1TP_142036_19910403_20170127_01_T1	30 m	5.00	Flow line delineation
Landsat TM	LT05_L1TP_142036_19930830_20170117_01_T1	30 m	0.00	Flow line delineation
Landsat TM	LT05_L1TP_142036_19980524_20161224_01_T1	30 m	0.00	Outline delineation
Landsat ETM+	LE07_L1TP_142036_20000724_20170210_01_T1	30 m	30.00	Outline delineation and flow line delineation
Landsat ETM+	LE07_L1TP_142036_20001028_20170209_01_T1	30 m	3.00	Outline delineation
Landsat TM	LT05_L1TP_142036_20060919_20161119_01_T1	30 m	1.00	Flow line delineation
Landsat ETM+	LE07_L1TP_142036_20120911_20161129_01_T1	15 m	16.00	Outline delineation
Landsat ETM+	LE07_L1TP_142036_20120927_20161128_01_T1	15 m	0.00	Outline delineation and flow line delineation
Landsat OLI	LC08_L1TP_142036_20150912_20170404_01_T1	15 m	0.28	Outline delineation, flow line delineation, and velocity measurement
Landsat ETM+	LE07_L1TP_142036_19991111_20170216_01_T1	15 m	7.00	Velocity measurement
Landsat ETM+	LE07_L1TP_142036_20011116_20170202_01_T1	15 m	1.00	Velocity measurement
Landsat ETM+	LE07_L1TP_142036_20101109_20161212_01_T1	15 m	2.00	Velocity measurement
Landsat ETM+	LE07_L1TP_142036_20101125_20161211_01_T1	15 m	1.00	Velocity measurement
Landsat ETM+	LE07_L1TP_142036_20111230_20161204_01_T1	15 m	0.00	Velocity measurement
Landsat ETM+	LE07_L1TP_142036_20120131_20161204_01_T1	15 m	2.00	Velocity measurement
Landsat OLI	LC08_L1TP_142036_20131227_20170427_01_T1	15 m	1.69	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20140723_20180524_01_T1	15 m	1.28	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20140909_20170419_01_T1	15 m	33.79	Velocity measurement
Landsat OLI	LC08_L1TP_142036_20151014_20170403_01_T1	15 m	8.63	Velocity measurement
Landsat OLI	LC08_L1TP_142036_20160914_20170321_01_T1	15 m	11.82	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20171003_20171014_01_T1	15 m	0.40	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20180920_20180928_01_T1	15 m	10.83	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20190907_20190917_01_T1	15 m	1.05	Flow line delineation and velocity measurement
Landsat OLI	LC08_L1TP_142036_20200909_20200918_01_T1	15 m	0.89	Flow line delineation and velocity measurement

.

2.2. Methods

2.2.1. Glacier Outline Delineation and Glacier Length

For the glacier outline in 2000, 1 Landsat TM image in 1998 and 2 Landsat ETM+ images in 2000 were selected. For the glacier outline in 2012, 2 Landsat ETM+ images were selected. For the glacier outline in 2015, 1 Landsat OLI image was selected. We manually digitized glacier outlines with the false-color composition (RGB by bands 543 for TM and ETM+ images and RGB by bands 654 for OLI images). To improve the resolution, multispectral Landsat ETM+ and OLI images were fused with the corresponding panchromatic band before digitization. The geometric union of glacier outlines in 2000 and 2012 was treated as one glacier inventory (GI_ST) to represent the glacier extent for DEM differences between SRTM and TDM90. The geometric union of glacier outlines in 2000 and 2015 was treated as one glacier inventory (GI_SH) to represent the glacier extent for DEM differences between SRTM and HMA. For glacier G085885E34389N, we manually digitized its main flow line from the 1971 topographic map (before digitization, we co-registered this topographic map to the 2006 Landsat TM image, and the co-registration error was 7.3 m). The main flow-line digitization followed these rules: start at the highest point, end at the lowest point, down along the center of the glacier [48]. Main flow lines for other dates were extended or shortened at the terminus according to the corresponding Landsat images. The length of the main flow line was taken as the length of glacier G085885E34389N.

2.2.2. Datum Adjustment

To compare elevations in different DEMs, they should use the same reference surface. We adjusted the ellipsoidal heights of SRTM X DEM, HMA, and TDM90 to the same orthometric heights as SRTM C DEM [49]

$$\mathrm{H}=\mathrm{h}-\mathrm{N},$$

(1)

where H is the orthometric height, h is the ellipsoidal height, and N is the geoid undulation between the EGM96 geoid and the WGS84 ellipsoid. The EGM96 geoid model from the National Geospatial-Intelligence Agency (NGA, https://earth-info.nga.mil, accessed on 4 February 2021) was used to implement this adjustment. After this adjustment, all DEMs were projected to UTM 45N and resampled to the resolution of 30 m by applying a bilinear interpolation.

2.2.3. DEM Co-Registration

Horizontal and vertical offsets exist widely among DEMs from different sensors, creating misalignment in DEMs. These horizontal and vertical offsets can be corrected using the triangular relationship of elevation differences (dh) with aspects and slopes on the stable area (without glaciers and lakes) [50]. Lakes in our study area were digitized from the Landsat OLI image in 2015, considering that the lakes in the Inner Tibetan Plateau are expanding [51]. TDM90 and HMA were co-registered to SRTM C DEM by applying the respective triangular relationships. During the co-registration, pixels with slopes >30° and <5° as well as elevation differences >100 m were taken as anomalies and discarded. We stopped the iterated co-registration when the improvement in the standard deviation (STD) of elevation differences was less than 10%. Elevation difference maps before and after co-registrations demonstrated that TDM90 and HMA were co-registered well to SRTM C DEM (Figure A3). To test the internal consistency of co-registration corrections, we also co-registered HMA to TDM90. Iteration times and displacements are shown in

Table 1. Parameters of DEM co-registrations.

Vector	Iteration Times	Displacements
		X (m)	Y (m)	Z (m)
HMA→SRTM	2	10.78	10.74	1.63
TDM90→SRTM	2	16.36	6.65	2.08
HMA→TDM90	1	−4.44	10.58	−0.74

. Given the good triangulation among the correction vectors of these 3 DEMs, our co-registration correction was found to have good internal consistency.

2.2.4. Surface Trend, Curvature, and Penetration Corrections

We fitted a second-order surface trend on the stable area in the elevation difference map between HMA and SRTM C DEM (Figure A4), potentially related to the tilt of the sensor such as that detected in the elevation differences between SRTM C DEM and KH-4/KH-9-derived DEMs [52,53,54]. It was eliminated by subtracting this trend.

High-resolution DEMs provide a better representation of terrain than low-resolution ones. Elevation differences from DEMs with different original resolutions can be determined using the maximal curvature [55]. Such errors can be compensated by the relationship between the elevation difference and the maximal curvature established over the stable area. The maximal curvature was computed on the finer DEM using the open-source software, the Geographic Resources Analysis Support System (GRASS GIS) (https://grass.osgeo.org/, accessed on 04 February 2021). We adopted this correction to the difference map between TDM90 and SRTM C DEM and it worked well (Figure A5). However, this correction did not obviously improve the difference map of HMA and SRTM C DEM (Figure A5), which may be attributed to the voids in HMA. Therefore, we dismissed this correction for the difference map of HMA and SRTM C DEM.

SRTM C DEM penetrates further into snow and ice than SRTM X DEM due to the different frequencies (5.7 GHz for C-band and 9.7 GHz for X-band) [55]. According to the SRTM configuration and data processing, the Attitude and Orbit Determination Avionics (AODA) was used to provide information regarding baseline length, attitude, and position for compensating the mast oscillation error [56,57]. However, residual motion effects remain after this oscillation correction, especially in SRTM X DEM, since the error can be slightly reduced in SRTM C DEM by the multiple coverages of ascending and descending paths [58]. This error only occurs in the along-track direction (personal communication with Michael Eineder). However, we did not detect this difference in ZK. Hence, we used the differences on glaciers as differences in penetration depth between SRTM X DEM and SRTM C DEM (Figure A6). Elevation differences between SRTM C DEM and SRTM X DEM were used to adjust the C-band to the X-band elevation values. A polynomial relationship built from the average elevation difference of 50 m elevation bins and the corresponding altitude on glaciers was applied to this correction (Figure A7). After this correction, the corrected SRTM C DEM approximately represents the glacier surface of later 1999. This approximation was taken into consideration when assessing the uncertainty of mass balance since X-band data also penetrate into snow and ice.

2.2.5. Glacier Mass Balance Estimation

We found erroneous elevation differences on the glacier between DEMs, which were attributed to the artificial elevation, the steep slope, or the shade. We segmented GI_ST and GI_SH into individual glaciers by ice divides derived from the second Chinese glacier inventory. Then, we adopted a sigmoid function to filter out elevation difference outliers, glacier by glacier, considering both the glacier elevation change always showing a non-linear profile [59] and the preservation of the inverse elevation change pattern of surge [60]

$$\mathrm{w}=(Ele{v}_{max}-Ele{v}_{min})/(Ele{v}_{gla}-Ele{v}_{min}),$$

(2)

$$r_{STD}=5-5\tanh (2\pi -5w),$$

(3)

$$\Delta h_{max}=r_{STD}\times ST{D}_{gla},$$

(4)

where $\mathrm{w}$ is the normalized elevation calculated from the maximal elevation ($Ele{v}_{max}$), the minimum elevation ($Ele{v}_{min}$), and the glacier elevation ($Ele{v}_{gla}$); $r_{STD}$ is the weight calculated from the sigmoid function; $\Delta h_{max}$ is the limit of elevation change and is the product of $r_{STD}$ and $ST{D}_{gla}$; the latter is the STD of the elevation differences on glaciers. Then, we needed to fill in the missing pixels. We selected 5800 m as the average ELA since 2000, considering both ELAs in the first Chinese glacier inventory and the average elevation difference in the difference map between HMA and SRTM C DEM (Figure A8). We allocated 0 to missing pixels at the accumulation zone (above 5800 m) with a slope >30° to minimize the slope-oriented error [52]. Other missing pixels were filled by the local (individual glacier) mean of the 100 m bands, as this method performs better compared with other interpolation methods [61]. The filled elevation differences were transferred to mass balance by multiplication with the snow/ice density of 850 ± 60 kg m⁻³ [62].

2.2.6. Glacier Velocity Extraction

The annual surface velocities of glacier G085885E34389N were measured from the panchromatic bands of 6 Landsat ETM+ images from 1999 to 2012 and 8 Landsat OLI images from 2013 to 2020 using the COSI-Corr software package [63,64], which has been widely used for displacement extraction of glaciers [65,66]. We built a subset covering glaciers in ZK to minimize the processing time. The Scan Line Corrector (SLC) of Landsat ETM+ failed on 31 May 2003. Data in SLC-off images lost ~22%. A weighted linear regression (WLR) method [67] was adopted to recover the missing pixels. Images on 9 November 2010 and 31 January 2012 were used to repair images on 25 November 2010 and 30 December 2011, respectively. Horizontal displacements in the North/South direction ($D_{NS}$) and the East/West direction ($D_{EW}$) were provided by the frequency correlator. The following parameters were determined after many tests: the initial window size was 64 × 64 pixels and the final window size was 16 × 16 in pixels, the step of the sliding window was 8 pixels, the robustness iteration was 4, and 0.95 was set as the mask threshold to obtain reliable results with a final resolution of 120 m. Glacier surface displacement (D) was calculated from two orthogonal components:

$$D=\sqrt{D_{NS}^2+D_{EW}^2}.$$

(5)

Then, the annual velocity (V, in m a⁻¹) was calculated considering the time interval (d, in days) between image pairs:

$$\mathrm{V}=D\times 365/\mathrm{d}.$$

(6)

Then, we smoothed the annual velocity using the Gaussian Blur function with a kernel size of 3 × 3 (360 × 360 m) and a standard deviation of 3. The annual velocities with values >600 m a⁻¹ were discarded. As a consequence, eight annual velocity results were acquired.

Velocities from the chosen GoLIVE NetCDF files with values >5 m d⁻¹ were taken as outliers and discarded. Then, they were converted to annual velocities (unit: m a⁻¹) by multiplication by 365.

The velocity profiles of glacier G085885E34389N were extracted along the nominal flow line shown in Figure 1b to highlight the evolution of the velocity.

2.2.7. Uncertainty Assessment

The delineation errors of glacier inventories in 2000, 2012, and 2015 can be calculated from glacier outlines and image resolutions [39]

$$se(S_I)=\mathrm{N}\times {\lambda }^2/2,$$

(7)

where $se(S_I)$ is the area of the uncertainty of the individual glacier inventory derived from Landsat images, N is the count of pixels along the glacier outlines, and $\lambda$ is the resolution of images (30 m for Landsat TM images and 15 m for Landsat ETM+ and OLI images). Area uncertainties were ±2.43%, ±1.32%, and ±1.44% of the glacier inventory in 2000, 2012, and 2015, respectively. Area uncertainties of GI_ST and GI_SH were calculated by the propagation of random errors

$$se(\stackrel{·}{S})=\sqrt{se{(S_i)}^2+se{(S_j)}^2},$$

(8)

where $se(S_i)$ and $se(S_j)$ are area uncertainties of different individual glacier inventories. For example, the area uncertainty of GI_ST was calculated from the uncertainties of glacier inventories in 2000 and 2012.

The uncertainties in the glacier length were due to the digitization at the terminus as the starting point with a high altitude was unchanged while the endpoint with a low altitude (at the terminus) always changed. We assumed errors were ±27 m for the glacier length from the 1971 topographic map [39] and half of the resolution for those from Landsat images (±30 m for MSS images, ±15 m for TM images, and ±7.5 m for ETM+ and OLI images).

The uncertainty of glacier mass balance from 1999 to 2015 was estimated in line with [68]

$$se(\stackrel{·}{M})=\sqrt{{(se(h_{pen})\times \rho /\Delta T)}^2+M^2\times ({(se(\rho )/\rho )}^2+{(se(\stackrel{·}{S})/S)}^2+{(se(\stackrel{·}{h})/\Delta h)}^2)},$$

(9)

including the penetration uncertainty ($se(h_{pen})$), the glacier density $\rho$, the time span ($\Delta T$), the calculated mass balance ($M$), the density uncertainty ($se(\rho )$), the glacier area uncertainty ($se(\stackrel{·}{S})$), the glacier area ($S$), the elevation change uncertainty ($se(\stackrel{·}{h})$), and the average elevation change ($\Delta h$). We assumed $se(h_{pen})$ was 1 m, in line with that in WKM [69], which is larger than the penetration depth of 0.61 ± 0.04 m in Puruogangri [70]. The density uncertainty was ±60 kg m⁻³ according to our density scenario. $se(\stackrel{·}{S})$ can be calculated from Equations (7) and (8). The elevation change uncertainty ($se(\stackrel{·}{h})$) was estimated over the stable terrain [71] as

$$se(\stackrel{·}{h})=E_{\mathsf{\Delta }{h}_i}/\sqrt{N_{eff}},$$

(10)

where $E_{\Delta h_i}$ equals to the STD of the mean elevation changes of the corresponding altitude bands and $N_{eff}$ is the effective count of independent elevation difference pixels

$$N_{eff}=N_{tot}\times P/2d,$$

(11)

where $N_{tot}$ is the total count of elevation difference pixels, $P$ is the DEM resolution (30 m), and $d$ is the distance of the spatial autocorrelation decided by Moran’s I autocorrelation index. In this study, d was 480 m for the elevation difference map between TDM90 and SRTM C DEM and 380 m for the elevation difference map between HMA and SRTM C DEM. For the mass balance from 1999 to December 2011 (from SRTM DEM and TDM90), there was another time-related uncertainty. We calculated the mass balance under different scenarios: TDM was acquired solely in 2011 or 2012. Thus, the calculation of mass balance from 1999 to December 2011 can be:

$$se(\stackrel{·}{M})=\sqrt{{(se(h_{pen})\times \rho /\Delta T)}^2+M^2\times ({(se(\rho )/\rho )}^2+{(se(\stackrel{·}{S})/S)}^2+{(se(\stackrel{·}{h})/\Delta h)}^2)+{M_t}^2}.$$

(12)

The calculated time-related mass balance uncertainty ($M_t$) was ±0.002 m w.e.a⁻¹, which means that the time-related uncertainty had a very limited influence on the final mass balance uncertainty.

Uncertainties in velocities were estimated as the mean velocity on the stable terrain. The mean uncertainty of velocities derived from COSI-Corr is 19.84 m a⁻¹. Velocity uncertainties were 56.67, 40.92, 16.64, and 26.22 m a⁻¹ for GoLIVE data with 16-, 32-, 80-, and 96-day intervals, respectively.

3. Results

3.1. Glacier Elevation Change and Mass Balance

The filled elevation difference maps are shown in Figure 2. Elevations decreased intensively at the lower part of glaciers, while the dominant thickening occurred at the high-altitude accumulation zone from 1999 to 2015. Under the density scenario of 850 ± 60 kg m⁻³ and our assessment of uncertainties, glacier mass balance was estimated to be +0.02 ± 0.07 m w.e.a⁻¹ from 1999 to 2011/12 and +0.03 ± 0.06 m w.e.a⁻¹ from 1999 to 2015. The eastward-flowing glacier, G085885E34389N, showed an inverse elevation change pattern: the reservoir area thickened and the receiving area thinned in the difference map between HMA and SRTM, as well as the difference map between ASTER DEMs and SRTM (Figure 2b–d), which can be a signal for a surge. Besides, the lowest part of glacier G085885E34389N was thinning until 2015, which represented that the surging mass had not reached the margin.

3.2. Other Surging Characteristics

There were crevasses on the reservoir area of glacier G085885E34389N in the Landsat OLI image in July 2014 (Figure 3a). Through the propagation of the surge, these crevasses moved down the glacier and the shear margin appeared at the southern edge in 2015 (Figure 3b). The terminus changed to advancing since 2016 and a periglacial lake was covered by the glacier in September 2020 (Figure 4a). There was an advance in the 1970s according to the glacier length change (Figure 4b; glacier outlines from 1971 to 2013 can be seen in Figure A9). The glacier has retreated since the 1980s, accelerating gradually to >20 m a⁻¹ after 2000. The retreat slowed down during 2014–2016. The glacial terminus transitioned to advancing, which lasted until September 2020. The terminus advanced ~510 m from September 2016 to September 2020. The increase of glacier length from 1971 to 1977 might be an earlier surging event.

3.3. Velocity of Glacier G085885E34389N during the Surge

According to the annual velocities derived from the COSI-Corr package (Figure 5a), we identified an apparent downward motion of a surging front. The velocities of glacier G085885E34389N increased from 2013 to 2015 and decreased after 2017. However, it was hard to determine when the motion reached the maximum. GoLIVE data with the 96-day interval showed that the velocity in the winter of 2015–2016 was higher than in the winter of 2016–2017 (Figure 5b). The slow-down between 2016 and 2018 was faster than that during 2018 and 2020. GoLIVE data with the 32-day interval showed that the glacier reached the fastest flow in October–November 2015 (585.77 m a⁻¹ at ~8000 m along the nominal flow line, Figure 5c). A sharp acceleration occurred during the summer of 2015 with an apparent downward motion of the surging front. Velocities from 2018 to 2019 (Figure 5d) showed a similar pattern to those from 2015 to 2016: the highest velocity occurred in October–November, and there was a deceleration in the winter.

4. Discussion

4.1. DEM Accuracy

As mentioned in [50] and done in [69], we downloaded Ice, Cloud, and land Elevation Satellite (ICESat) data from NSIDC for the assessment of DEM accuracy. We adopted the elevation correction in [49] to adjust the ICESat heights relative to the EGM96 geoid and the co-registration in [50] before the assessment. During the assessment, points with elevation differences of ≥100 m relative to DEMs and a slope >15° (to ensure the ICESat had an accuracy of ≤1 m [72]) were discarded. The average and root mean square error (RMSE) of the elevation difference over the stable area were 0.42 and 2.53 m, 0.15 and 2.24 m, and 0.02 and 0.91 m between three DEMs (SRTM C DEM, TDM90, and HMA) and ICESat data, respectively.

4.2. Mass Balance Comparison

The two mass balance results in this paper are consistent (+0.02 ± 0.07 m w.e.a⁻¹ from 1999 to December 2011 and +0.03 ± 0.06 m w.e.a⁻¹ from 1999 to 2015); hence, it is robust to some extent in that the mean annual mass balance in ZK was slightly positive/balanced from 1999 to 2015. The ICESat-derived mass balance in the ZK and Songzhi Peak region was +0.37 ± 0.25 m w.e.a⁻¹ during 2003–2009 [41]. This discrepancy can be attributed to the gathering of ICESat altimetry points on the accumulation zone (Figure A10), which led to an erroneous result. In the early 2000s, some glaciers on the Karakoram were observed to be thickening [73]; this phenomenon came to be known as the “Karakoram anomaly”. Later, regional mass balance studies employing ICESat laser altimetry points confirmed this mass gain on the Karakoram, while more significant glacier mass gain occurred on the Central and Northern part of the TP, which suggested the center of the glacier mass anomaly may lie in the inner part of the TP [41,74]. With the enrichment of remote sensing data, systematic glacier mass balance assessments from sequential DEMs (TSX/TDX-derived DEMs, ASTER DEMs, and HMA DEMs) of the TP demonstrated that the most significant mass gain area among the TP was a strip from the eastern WKM to the upper Tarim Basin rather than the Karakoram, while glaciers in ZK were just nearly balanced from 2000 onward [42,43,75]. Considering the limited coverage of ICESat points on glaciers on the TP, the erroneous mass balance result caused by the gathering of ICESat points on ZK, and the consistent mass balance results of ZK in this paper, we suggest that glaciers in ZK have had a slightly positive/nearly balanced condition since 2000 and ZK may be the southern edge of this glacier anomaly on the TP.

4.3. Evolution of the Surge

The normal elevation change pattern in the difference map between TDM90 and the SRTM C DEM and the elevated low-altitude receiving area of glacier G085885E34389N in the elevation difference maps between two ASTER DEMs and SRTM C DEM (Figure 2) indicate that the surge occurred between October 2012 and May 2014. In terms of the definition of the glacier surge, however, rapid and abnormal velocity should be the principal characteristic of a glacier surge. Therefore, it is vital to determine when the velocity was normal and abnormal. According to the flow line change, the terminus retreat rate was relatively consistent during the 1990s and 2000s. The glacier velocity in this period was representative of normal/quiescent velocity. Unfortunately, we do not have Landsat OLI images in ZK produced before March 2013. The scan line corrector (SLC) failure caused massive data loss, presented as repeating along-scan stripes in Landsat ETM+ images since May 2003 [76]. The low 30-m resolution of Landsat TM images and thus the lack of surface characteristics make them unsuitable for velocity extraction. Two Landsat ETM+ images on 11 November 1999 and 16 November 2001 were used for measuring normal/quiescent velocity. Velocities from November 2010 to December 2011 were more akin to velocities from 1999 to 2001 (though partly affected by artificial pixel recovering), and velocities in the winter of 2013–2014 were higher than the quiescent velocity (Figure 6). The maximal velocity at 6000–9000 m along the nominal flow line in early February 2014 was over 200 m a⁻¹, while the velocity in 1999–2001 was under 20 m a⁻¹. Considering that the normal velocity should be slowest annually during winter, abnormal (over 10 times) velocity should be a mark of a surge. Moreover, the winter acceleration from December 2013 to February 2014 was an abnormal phenomenon since there were decelerations in the winters of 2015–2016 and 2018–2019 (Figure 5). This abnormal winter acceleration can also be related to a surge. The surge after October 2012 may have already started before January 2014. The maximal velocity in October–November 2015 indicated that the acceleration of the surge lasted roughly two years. The terminus started to advance in 2016, which was approximately three years later than the initiation and, notably, roughly one year later than the maximal velocity of this surge. The ongoing terminus advance and the high winter velocity of 2019–2020 compared with that in the winter of 2013–2014 indicated that the deceleration of the surge had been longer than five years. As a consequence, the active phase of this surge had been over seven years. The terminus advance in the 1970s (Figure 4b) indicated there may be a previous surge event, and the recurrent duration may be longer than 40 years.

4.4. Mechanism of the Surge

The long active phase of over 7 years, the shorter acceleration phase of ~2 years compared with the deceleration phase of over 5 years, and the possible long surge duration of over 40 years are more similar to the Svalbard-type surge rather than the Alaskan-type surge. The relatively low velocity of 587 m a⁻¹ and the slow advance of terminus are quite different from the Alaskan-type surging glacier such as the Variegated Glacier [7] and the Sortebræ Glacier [77]. However, the winter acceleration during the acceleration phase from December 2013 to February 2014 was akin to the Alaskan-type surge [7,10,11,77]. The maximal velocity of a calendar year, which occurred in October–November in ZK, was similar to the West Kunlun Glacier and the N2 Glacier in WKM, which showed peak velocities from the late fall to winter [17]. Yasuda and Furuya attributed this peak velocity of these two glaciers in WKM to the lubrication of remanent summer meltwater and the reduction of macroscopic porosity during the late fall to winter, which is akin to the hydrologically controlled Alaskan-type surge. The same processes may also occur in ZK. Therefore, we cannot classify this surging glacier as either of the two well-known types of surges.

Some other glaciers among the Central and Northern TP have shown surging characteristics in the Muztag [32], Western Kunlun Mountains (WKM) [5,36], Puruogangri [33], Geladaindong [35], and Xinqingfeng and Malan ice caps [34] since 2000. The surging glacier in Muztag belongs to the Eastern Kunlun Shan region [43], which showed a balanced mass condition of −0.07 ± 0.07 m w.e.a⁻¹ from 2000 to 2018. Glaciers in WKM have gained mass since the 1970s [69]. The surging glacier in Puruogangri [33] showed a slight mass gain of +0.079 ± 0.015 m w.e.a⁻¹ (+0.016 ± 0.023 m w.e.a⁻¹ using a different method) from 2000 to 2012 [78]. Surging glaciers in Geladaindong experienced a relatively balanced glacier mass change (−0.22 ± 0.08 to +0.09 ± 0.03 m w.e.a⁻¹) from 2000 to 2012 [79]. A similar mass condition (−0.29 ± 0.18 to 0.10 ± 0.16 m w.e.a⁻¹) was reported for surging glaciers in the Xinqingfeng and Malan ice caps from 1999 to December 2011 [34]. These mass balance results indicate considerable thickening in the reservoir part. We checked the area, length, and slope in Randolph Glacier Inventory 6.0 (RGI 6.0, available from http://www.glims.org/RGI/, accessed on 4 February 2021) at six surging regions among the Central and Northern TP. The surging glaciers have geometric similarity in terms of their large size and flat surface (

Table 2. Averaged geometric information of glaciers in six surging regions.

Region	Category	Area (km2)	Slope (°)	Length (km)
	Surge	35.6	7.9	13.5
	Non-surge	3.4	22.3	2.2
Western Kunlun Mountains (WKM)	Surge	96.0	9.2	19.1
	Non-surge	4.8	22.8	2.4
Puruogangri	Surge	14.1	12.0	7.6
	Non-surge	6.5	18.0	3.2
Geladandong	Surge	26.1	11.5	10.2
	Non-surge	5.0	18.3	2.9
Xinqingfeng and Malan	Surge	57.2	8.2	14.5
	Non-surge	3.7	18.3	2.4
Muztag	Surge	96.6	14.3	13.2
	Non-surge	3.3	20.7	2.0

), which are coincident with the observations: longer and flatter glaciers tend to surge [22]. A dynamic—geometric relationship also affirms that a surging glacier has to exceed a “critical length”, which means longer glaciers are more likely to surge [16]. According to the enthalpy balance theory [19,80], the thickening in the reservoir part leads to the increase in enthalpy, while the large size and gentle slope impede the loss of basal water, which may jointly trigger a surge. Hence, the relatively balanced mass condition, large area, and shallow slope may act as controls on glacier surges among the Central and Northern TP.

Many other elements can control surges. To understand the surge thoroughly and identify a glacier surge in a timely manner, regular in situ and remote-sensing-based monitoring of climate, glacier motion, thermal regime, drain system, and basal condition are required. A comprehensive interpretation of the dynamic mechanism of surging glaciers on the Central and Northern TP is still being developed.

5. Conclusions

Our study revealed that glaciers in ZK were under a slightly positive/nearly balanced mass condition of +0.03 ± 0.06 m w.e.a⁻¹ from 1999 to 2015 (+0.02 ± 0.07 m w.e.a⁻¹ from 1999 to December 2011). This glacier mass condition and other studies of glacier mass balance on the Central and Northern TP may suggest ZK is located at the southern edge of the glacier anomaly over the TP. We found an eastward-flowing glacier (GLIMS ID: G085885E34389N) that showed surging signals, such as inverse elevation change pattern, crevasses, shear margin, and snout advance. Sequential DEMs and velocity data indicate that the surge started between October 2012 and January 2014. There was a possible surging event in the 1970s on the same glacier. This surging glacier found on the Central TP, as well as other published surging glaciers among the Central and Northern TP, may extend the territory of surging glaciers on the TP. The relatively positive mass budget, large size, and shallow slope, however, may control the occurrence of surges in this area. The explicit surging mechanism accounting for surges on the Central and Northern TP needs more exploration, as the lack of in situ observations and information on the glacial-thermal regime and the basal condition hindered our interpretation.