Morphometric Controls on Glacier Mass Balance of the Puruogangri Ice Field, Central Tibetan Plateau

Abstract

Evaluating the impacts of climatic changes and morphometric features on glacier mass balance is crucial to providing insight into glacier changes and their effects on regional water resources and ecosystems. Here, we presented an evaluation of morphometric effects on the glacier mass balances of the Puruogangri ice field (PIF) on the Tibetan Plateau. A clear spatial variability of glacier mass balances, ranging from −0.035 to +0.019 m·w.e.·year⁻¹, was estimated by comparing the TanDEM-X DEM (2012) with the SRTM-X DEM (2000). In general, the observed glacier mass changes were consistent with our fieldwork investigations. Furthermore, by applying the method of linear regression analysis, we found that the mass changes of individual glaciers on the PIF were mainly dominated by the mean altitude (R = 0.84, p < 0.001), however, they were statistically independent of glacier size, aspect, and surface velocity. At a local scale (grid size of 10 × 10 pixels), apart from the factor of altitude, surface velocity was correlated with glacier mass change.

1. Introduction

Mountain glaciers on the Tibetan Plateau (TP) and its surroundings play a key role in the water and food security of people living in eastern and southern Asia because these glaciers are located at the headwaters of many prominent Asian rivers [1,2]. Therefore, it is necessary to identify the factors controlling glacier mass balance and provide insight into glacier changes and their effects on regional water resources and ecosystems. Basically, the mass balance of mountain glaciers is controlled by regional climate variability and local morphometric factors [3]. In recent decades, our knowledge of mountain glacier responses to climate changes has been substantially improved with advances in regional climate models (RCMs) and climatic observations [4,5,6,7]. Moreover, the spatially distributed glacier mass balances, which were mainly derived with geodetic method [8,9,10,11,12,13] and glacier surface energy and mass balance model [14,15,16,17,18], provide opportunities for investigating the effect of morphometric factors on glacier mass changes.

In general, previous research suggested that there are several effective morphometric factors of glacier mass balances on the TP and its surroundings, such as topographic characteristics (i.e., glacier size, altitude, surface slope, curvature and aspect) [19,20,21], ice surface dynamics [22,23], and debris covers [24,25]. Racoviteanu et al. [21] reported that glacier changes in the Kanchenjunga-Sikkim area (eastern Himalaya) were primarily controlled by glacier sizes and altitudes. Berthier et al. [19] found that glacier mass balance in the Himachal Pradesh (western Himalaya) was mainly affected by glacier altitudes and not correlated with size and aspect. In terms of glacier dynamics, Bolch et al. [22] and Gardelle et al. [23] agreed that low glacier surface velocity may result in a clear mass loss; however, a relatively stable front was observed over the stagnant tongue regions of glaciers covered with heavy surface debris in the central Himalayas [24]. Furthermore, Gardelle et al. [23] found variable effects of debris cover on glacier mass balance over nine study sites distributed in the Pamir, Karakoram, and Himalayas, which might be dependent on the fraction and thickness of debris, as well as the glacier velocity [24]. Overall, the local-scale behavior of glacier mass changes is influenced by multiple morphometric factors and is not fully understood on the Tibetan Plateau and its surroundings [2], particularly on the central TP, where few investigations have been conducted.

In this study, we aim to evaluate the morphometric effects on glacier mass balances of the Puruogangri ice field (PIF) on the central TP, which is the largest modern ice field on the Tibetan Plateau. The glacier elevation changes and mass balances of the PIF between 2000 and 2012 were estimated by comparing the SRTM DEM (2000) with the TanDEM-X DEM (2012). Importantly, the observed glacier elevation changes were validated with our fieldwork investigations. We calculated glacier size from a well-delineated glacier outline derived from a Landsat TM image and extracted topographic terrain features (altitude, surface slope, curvature, and aspect) from the high-resolution TanDEM-X DEM. Ice flow velocity was derived according to the repeat-pass differential SAR interferometry method with a pair of ALOS PALSAR data. Furthermore, we employed linear regression analysis to evaluate the effects of topographic features and surface ice dynamics on mass balances of individual glaciers on the PIF. Additionally, the entire ice field was gridded with a size of 10 × 10 pixels to further estimate the morphometric effect on glacier mass changes at a more local scale.

2. Geographical Setting of the Study Area

The PIF is located on the northern central Tibetan Plateau and the western portion of the Tanggula Mountains, with central coordinates of 33°55′ N and 89°10′ E [26] (Figure 1). The glaciated area covers more than 420 km² and a volume of approximately 52 km³ [27]. The equilibrium-line altitude (ELA) of the PIF lay at 5748 m·a.s.l. between 2001 and 2012 [28]. The glacier system is composed of approximately 50 glaciers with variable lengths radiating from the center of the ice field and extending halfway down-valley to the piedmont terminus. Glacier melt-water contributes to three endorheic lakes, named Linggo Co, Dogai Coring, and Chibuzhang Co [29]. According to the lake catchments, the ice field is divided into three parts—i.e., the northern, eastern and western slopes.

The climate of the PIF is mainly affected by continental climatic processes, although it lies close to the boundary between the monsoon-dominated south and the westerly dominated north [30]. The records of the two ice cores drilled in 2000 (locations are shown in Figure 1) indicated an abrupt onset of the recent warming in the early 20th century [30]. Following the Little Ice Age, the northern slope glaciers demonstrated an unambiguous and abrupt retreat, whereas the terminus changes of glaciers on the eastern and western slopes were relatively more stable [29]. The geodetic measurements of surface elevation changes on the PIF exhibited a mean value of approximately −0.238 m·year⁻¹ from 1974 to 2000 [31] and −0.049 m·year⁻¹ between 2000 and 2012 [12]. Huintjes et al. estimated glacier elevation change of −0.58 m between 2000 and 2011 by using the ‘COupled Snowpack and Ice surface energy and MAss balance model’ (COSIMA) [14]. The latest study revealed a much more rapid glacier thinning (−0.317 m·year⁻¹) for the period 2012–2016 than in the periods 1974–2000 and 2000–2012 [13]. This finding indicates an accelerating trend of glacier surface thinning over the PIF since the beginning of the 2010s.

3. Data and Methods

3.1. Glacier Mass Balance Estimation

3.1.1. TanDEM-X DEM and Glacier Mass Balance

We employed the DEM differencing method to calculate the elevation changes and mass balances of all glaciers over the PIF by comparing the TanDEM-X DEM (2012) and the SRTM-X DEM (2002). The TanDEM-X DEM of the entire ice field was derived from two sequential pairs of TerraSAR-X/TanDEM-X SAR imagery (their coverage is shown in Figure 1), which were acquired on 26 January 2012. The TanDEM-X bistatic SAR data were received as co-registered Single Look Complex (SLC) data in the CoSSC format, with 3.3 m and 2.5 m resolution pixel spacing in the azimuth and range directions, respectively. These original SLC data were processed according to the bistatic SAR interferometric method [13] to generate a glacier surface DEM of the PIF with a spatial resolution of 10 m. It is noted that absolute and relative height errors in the bistatic interferometric processing were corrected with independent height reference datasets derived from 1:50,000 topographic maps, ICESat laser altimeter data, and GPS measurements. The locations of these height reference data are shown in Figure 1.

The SRTM-X DEM only covers a large part of the ice field (~90%) because the X-band radar was operated to derive data along discrete 50 km wide swaths [32]. The remaining area of the SRTM-X DEM was filled by the 1 arc-second SRTM-C DEM (30 m spatial resolution) released in January 2015, and the penetration bias of the C-band and X-band radar wave for ice and snow was corrected following the method described by Gardelle et al. [23]. Furthermore, the 30 m resolution SRTM-X DEM was resampled using the bi-linear algorithm to 10 m resolution pixels to match with the TanDEM-X DEM.

Before the process of DEM differencing, the co-registration method of horizontal shift, which was suggested by Berthier et al. [19], was used to eliminate the influence of geometric errors induced by the horizontal difference between the two DEMs used [33]. In this proposed method [19], the shift values between the TanDEM-X DEM and SRTM-X DEM were determined by the minimum in the standard deviation of the two DEM elevation differences of the off-glacier regions. Once accurately co-registered, the DEMs were differenced to map elevation changes for individual glaciers over the PIF, in which the outlines of individual glaciers were extracted from Landsat TM images, as described in Section 3.2.

The resulting elevation changes were multiplied with the density of lost or gained glacial material between two periods of investigation to estimate the mass balances of all glaciers on the PIF. Because there are no measurements of the glacier ice density available in our study area, we employed an average glacial density of 850 kg·m⁻³. This value of glacial density recommended by Huss (2013) was estimated by applying an empirical firn-compaction model, and was proposed for geodetic observations with a time interval of more than five years [34].

3.1.2. Uncertainty Analysis

In this study, we applied a common method for accuracy assessment of elevation changes obtained by the DEM differencing [35] due to the lack of in situ observations available. This method is generally used to calculate the statistical error ($\sigma$) of the elevation differences over off-glacier regions according to the following equation:

$$\sigma =\sqrt{m^2+{(\frac{STDV}{\sqrt{N}})}^2},$$

(1)

where $m$ and $STDV$ are the mean value and standard deviation of elevation differences over off-glacier areas, respectively. N is the number of independent measurements and can be calculated with the Equation (2):

$$N=\frac{N_{total}\times PS}{2D},$$

(2)

where $N_{total}$ is the total number of measurements, PS is the pixel size of the DEMs used, and D is the distance of spatial autocorrelation. According to Bolch et al. [22], the decorrelation distance (D value) is generally equal to 20 pixels. Here, we utilized the 30 m resolution DEM and a decorrelation distance of 600 m.

Furthermore, we added an uncertainty of 7% when the volume change was converted into mass balance due to the lack of in situ measurements of the glacial density [36]. Thus, an error of ±60 kg·m⁻³ was used as the density conversion factor.

3.2. Glacier Boundary Delimitation

A well-delineated glacier outline is necessary to accurately estimate glacier mass balance and size. Glacier outlines can be downloaded from the Global Land Ice Measurements from Space (GLIMS) database of glacier inventory [37]. However, these boundaries do not always match our requirements of accuracy [23].

Satellite optical images—such as Landsat TM/ETM+, ASTER, and SPOT—have been demonstrated to be valuable data sources for delineating glacier outlines [23,24]. In this study, we applied the Normalized Difference Snow Index (NDSI), which has been used by Silverio and Jaquet [38], for a Landsat TM image acquired on 24 September 2007. Following the principle of the NDSI, glaciated areas show high values due to the snow and ice coverage, whereas areas of no snow coverage in the surrounding regions show low values. According to the histogram of NDSI images, a threshold of 0.55 was used to delineate the glacier outlines. Moreover, the glacier boundaries were manually corrected for several glaciers to pursue high accuracy and precision. Finally, the sizes of all glaciers on the PIF were calculated from the well-delineated glacier outlines.

In order to validate the surging event of glacier 5Z213E0012, we detected changes in glacier terminus positions from 1976 to 2016, which were retrieved by one Landsat Multi-Spectral Scanner (MSS) image in 1976, two Landsat Thematic Mapper (TM) images in 1998 and 2010, two Landsat Enhanced Thematic Mapper Plus (ETM+) images in 2000 and 2002, and one Landsat Operational Land Imager (OLI) image in 2016.

3.3. Extraction of Topographic Parameters

The topographic terrain features (altitude, surface slope, curvature, and aspect) were extracted from the resulting TanDEM-X DEM with 10-m resolution, as described in Section 3.1. With regard to geometry and geography, the altitude is defined as the height above sea level, the slope is the feature used to describe the steepness of a terrain surface, curvature is the indicator of a surface’s roughness, and aspect is the bearing or azimuth of the slope direction. The altitude is equal to the elevation of the glacier surface DEM, whereas the other three topographic parameters require retrieval with GIS (geographic information system) software. In this study, we employed the spatial analyst toolbox in the ArcGIS software program (version 10.0) (http://www.esri.com/software/arcgis) to estimate the parameters of slope, curvature, and aspect. Moreover, mean values of these topographic parameters for individual glaciers over the PIF were estimated with the glacier outlines described in Section 3.2.

3.4. Ice Flow Velocity Determination

Speckle/feature tracking and repeat-pass differential SAR interferometry (D-InSAR) are two well-established methods for estimating glacier velocity [39,40]. Here, we applied the D-InSAR method to estimate the ice flow velocities of glaciers on the ice field using the DIFF module in the GAMMA software (http://www.gamma-rs.ch/). This method is prone to decorrelation due to ice melting and movement, and the L-band radar of the longer wavelength has a better capacity to preserve coherence over glacier surfaces than C-band and X-band radars [41]. Therefore, we used a pair of L-band ALOS PALSAR data over the PIF, acquired on 20 January 2009 and 7 March 2009, to minimize the decorrelation influence on the glacier velocity determination. Due to a lack of in situ observations of glacier surface velocity in this study area, accuracy assessment was based on the assumption that the ice-free areas are stable (no movement) [42,43]. The detailed calculations are similar with that of elevation changes, as described in Section 3.1.2. However, the $m$ and $STDV$ in the Equation (1) were replaced by the mean value and standard deviation of surface displacements over ice-free areas, respectively.

3.5. Statistical Analysis

The linear regression analysis was conducted to evaluate the dependence of the glacier mass balance over the PIF on the selected topographic factors (size, altitude, surface slope, curvature, and aspect) and the ice flow velocity. Specifically, we employed simple linear regression (SLR) analysis to calculate the correlation coefficient (R), which describes the linear relationship between glacier mass balance and single independent morphometric factors. In addition, multiple linear regression (MLR) analysis, as an extension of the SLR, was employed to examine the relationship between the glacier mass balance and two or more morphometric factors by fitting a linear equation to observed data. In this study, we used the SLR and MLR tools in the SPSS (Statistical Product and Service Solutions) version 19.0 software package.

3.6. Field Works

Since 2000, we carried out five field-working campaigns in order to study glacier dynamics and glacial moraines of the PIF, although it is very difficult to reach. During these in situ investigations, the surface morphologies in glacier termini regions were monitored with photographic records. Furthermore, we performed GPS surveying on the front of glacier 5Z611A0006 with RTK technique in October 2012 to provide a high-accuracy elevation reference for bistatic SAR interferometric processing. The reference station of RTK was installed on a rock outcrop and surveyed with the Unistrong E650 GPS receiver for more than five hours. Moreover, the records of base GPS receiver were post-processed using differential GPS (DGPS) technique with the Lhasa station of CMONOC (Crustal Movement Observation Network of China) to obtain a result of height accuracy up to $\pm 0.1\mathrm{m}$. Another GPS receiver of the same type with base receiver was simultaneously employed for roving at 45 points that are apart from the reference station less than 3 km. In particular, the error of relative range between the floating and reference station was required to be better than 1 cm, in order to keep the high elevation accuracy of these GPS points. The locations of GPS measurements are shown in Figure 1.

4. Results

4.1. Glacier Elevation Changes

Figure 2 illustrates detailed elevation changes over the entire PIF generated by comparing the SRTM-X/C DEM in February 2000 with the resulting TanDEM-X DEM in January 2012. The glacier-wide mean elevation difference between 2000 and 2012 is −0.27 ± 0.18 m (or −0.023 ± 0.015 m·year⁻¹), which compares to −0.51 ± 3.71 m (DEM differencing method) observed by Neckel et al. [12] in the same period. This discrepancy is presumably due to difference in the coverage of used TanDEM-X bistatic SAR data and the vertical accuracy of generated TanDEM-X DEMs. In this study, the entire PIF was covered by two pairs of TanDEM-X images, whereas only one pair of TanDEM-X data (covers about 90% of the ice field) was applied by Neckel et al. [12]. Moreover, our DEM result has a better vertical accuracy of 1.91 ± 0.76 m (see Liu et al. 2016) than that of 3.24 ± 1.11 m reported by Neckel et al. [12], because we have corrected the absolute and relative height errors in the bistatic interferometric processing by using the well-distributed height reference datasets (their locations are indicated in Figure 1).

Additionally, Huintjes et al. estimated glacier elevation change of −0.58 m between 2000 and 2011 with the ‘COupled Snowpack and Ice surface energy and MAss balance model’ (COSIMA) by using the meteorological dataset from the High Asia Refined analysis (HAR) [14]. This difference of −0.31 m between our observed results and modeled glacier elevation changes is possibly attributed to two reasons. Firstly, there are different time periods of investigation (2000–2012 and 2000–2011) in the two studies. Secondly, precipitation scaling factor determined at Zhadang glacier (30°28.57′ N, 90°38.71′ E) was used for the PIF (33°55′ N, 89°10′ E) by Huintjes et al. [14]. However, according to Yao et al. [2], precipitation variations at Zhadang glacier and the PIF are not similar due to different latitudes. Importantly, the modeled glacier elevation changes are highly sensitive to precipitation scaling factor, as exhibited in Table 3 of Huintjes et al. [14].

As shown in Figure 2, clear surface thinning was observed in most glacier tongue regions, with a maximum elevation decrease of up to approximately −40 m, whereas a slight thickening was detected in the broad and flat interior regions (mostly in accumulation areas) above 5800 m·a.s.l. In general, this spatial pattern of glacier elevation changes over the PIF was similar with that revealed by Neckel et al. [12]. Particularly, it is worth noting that the surface elevation changes varied for individual glaciers. Most glaciers experienced predominant surface lowering of approximately 20–40 m at glacier terminal regions, e.g., glacier No. 5Z611A0008 and No. 5Z513B0018; however, several glaciers (e.g., No. 5Z611A0006) exhibited no obvious changes in elevation, except for one thickening glacier (No. 5Z213E0012) in the eastern PIF. Probably due to surging activities, this anomalous glacier was thickened by up to approximately 80 m at the terminus, whereas negative elevation changes were observed in the upper ablation areas. Figure 3 highlights the detailed features of surface elevation changes and altitudes along the central flow line of these four selected glaciers.

The DEM-derived elevation changes at glacier tongues were generally consistent with our fieldwork investigations. For example, the contradictory behavior at the tongue regions of two typical glaciers (No. 5Z611A0006 and No. 5Z611A0008) could be interpreted by their photographic records from 2000 to 2012, as shown in Figure 4 and Figure 5. The relatively stable glacier 5Z611A0006, as evidenced by elevation change results (see Figure 2 and Figure 3a), experienced almost no geomorphological changes (Figure 4). However, for glacier 5Z611A0008, intensive surface melting at the terminal areas was found in our field campaign in 2011 (Figure 5). The melt-water and ice cracking occurring in the snout indicated that glacier 5Z611A0008 exhibited significant mass loss or instability.

4.2. Glacier Surface Velocity

Figure 6 shows the InSAR-derived ice flow velocities of the Puruogangri ice field in 2009. The glacier-wide mean surface velocity is 2.25 ± 0.61 m·year⁻¹, and the maximum velocity is up to about 17 m·year⁻¹. In general, our results are similar with the observed glacier surface velocities in 1996 (Figure 1 in Huintjes et al. [14]). For example, obvious surface movements were detected by two studies in the glacier tongue areas of glacier 5Z611A0003 and 5Z611A0006, and in the upper regions of glacier 5Z611A0008, 5Z611A0009, 5Z513B0007, and 5Z513B0014. However, the rates of surface velocities were significantly changed over glacier 5Z213E0010 and 5Z513B0010 between 1996 and 2009. In particular, for glacier 5Z213E0012, we found that ice flow velocities over upper regions were much higher than that of glacier tongue areas in 2009, whereas relatively faster movements were observed over glacier tongue areas in 1996 [14].

4.3. Spatial Variability of Glacier Mass Balances

A slightly negative glacier mass balance of −0.019 ± 0.014 m·w.e.·year⁻¹ was measured over the entire ice field from 2000 to 2012. We observed clear spatial variability of glacier mass balances ranging from −0.035 to +0.019 m·w.e.·year⁻¹ in the study area.

Table 1. Mass balances of glaciers with a size of greater than 1 km² on the PIF (except for the surge-type glacier) and glacier size and mean values of other morphometric factors.

WGMS ID	Mass Balance (m·w.e.·Year−1)	Size (km2)	Altitude (m)	Slope (Degree)	Curvature	Aspect (Degree)	Velocity (m·Year−1)
5Z611A0002	0.017	3.51	5849.5	17.7	0.004	210.0	1.36
5Z611A0003	0.000	6.60	5750.3	10.8	0.008	234.9	3.40
5Z611A0005	0.009	5.13	5878.3	16.2	0.034	217.5	2.06
5Z611A0006	0.004	19.86	5833.4	8.7	0.010	216.9	3.83
5Z611A0008	−0.021	46.28	5758.6	5.5	0.000	214.2	2.72
5Z611A0009	−0.006	11.14	5810.4	9.8	0.005	237.6	2.90
5Z611A0010	0.003	50.18	5770.5	8.5	0.001	207.0	2.40
5Z611A0013	0.005	3.99	5824.2	17.0	−0.040	225.0	1.00
5Z611A0014	0.018	1.34	5990.6	13.7	0.026	167.8	1.21
5Z611A0015	0.006	14.85	5831.0	9.1	0.008	189.6	2.55
5Z213E0001	−0.018	6.08	5668.3	12.3	−0.011	152.7	0.99
5Z213E0004	−0.019	1.86	5622.3	10.7	−0.008	199.3	2.03
5Z213E0008	−0.017	1.79	5645.0	14.1	−0.038	104.1	0.98
5Z213E0009	−0.035	7.86	5627.2	10.7	−0.023	130.5	1.78
5Z213E0010	−0.008	32.27	5787.3	12.4	−0.004	140.3	2.23
5Z513B0005	0.002	8.33	5841.6	15.5	−0.032	134.1	1.42
5Z513B0007	0.019	14.34	5917.1	19.9	0.016	119.7	1.85
5Z513B0009	0.013	4.38	5830.6	21.0	0.056	80.7	0.86
5Z513B0010	0.008	13.31	5773.8	12.2	0.003	126.2	4.50
5Z513B0011	−0.025	1.19	5552.6	15.1	−0.029	266.7	0.74
5Z513B0013	−0.012	2.34	5685.5	17.1	−0.051	231.1	0.70
5Z513B0014	−0.004	28.85	5791.0	15.7	−0.009	173.2	1.90
5Z513B0015	−0.006	1.62	5780.1	25.0	0.069	87.5	0.93
5Z513B0016	−0.012	4.31	5720.5	14.0	−0.006	109.7	1.45
5Z513B0022	−0.013	6.17	5665.9	13.5	−0.019	210.0	0.76
5Z513B0020	−0.002	5.23	5767.3	18.5	−0.001	98.3	1.17
5Z513B0018	−0.007	53.58	5775.6	9.4	−0.007	145.9	1.68
5Z513C0001	−0.011	10.42	5657.3	10.6	−0.008	194.4	1.37
5Z513C0002	0.004	3.45	5701.3	17.4	0.034	266.8	0.52
5Z513B0024	−0.003	32.34	5689.2	16.1	−0.011	165.8	1.42

lists the mass balances, glacier sizes, and mean values of other morphometric factors for 30 glaciers on the PIF with a size of greater than 1 km² (except for the surge-type glacier). A total of 18 glaciers experienced mass loss, but mass gain occurred in 12 glaciers. Furthermore, a heterogeneous pattern of glacier mass balance was observed within three sub-regions corresponding to three topographic units of the northern, western and eastern parts in the PIF (see Figure 7). On the western slope, there were three glaciers with negative mass balance and seven glaciers with positive mass balance. Similarly, 10 mass-loss glaciers and 5 mass-gain glaciers were observed on the northern slope. However, on the eastern slope, the behavior of all five glaciers was relatively consistent, demonstrating distinct mass loss.

4.4. The Effect of Morphometric Factors on Glacier Mass Changes

The results of both SLR and MLR analysis revealed variable morphometric effects on the mass balance of individual glaciers on the PIF (see Figure 8 and

Table 2. Summary of the MLR analysis in the SPSS software package (T-value (t) and p-value are two parameters obtained from the t-test).

Variable	Parameter Estimate	Standard Error	t	p-Value
Size	0.002	0.010	0.174	0.863
Altitude	0.010	0.002	6.538	<0.001
Slope	0.118	0.049	2.393	0.025
Curvature	1.661	5.743	0.289	0.775
Aspect	0.005	0.003	1.822	0.082
Velocity	0.215	0.181	1.191	0.246

). In general, the relationships between glacier mass balance and the morphometric controls modeled by SLR analysis are consistent with those revealed by MLR analysis. The results showed a strong positive correlation (R = 0.84 > 0.75 in Figure 8f) between glacier mass balance and the mean altitude at the 99% confidence interval (p < 0.001 in Table 2). Glacier surface slope might be slightly correlated with glacier mass change (R = 0.35 > 0.25) in Figure 8d at the 95% confidence interval (p = 0.025 < 0.05). Furthermore, the SLR analysis suggested that glacier size, aspect, and ice flow velocity were not statistically significant controls on glacier mass balance (see Figure 8a–c). These findings were generally consistent with the p-values of the above three factors (p > 0.05, 95% confidence interval) in Table 2. However, the possibly positive correlation between curvature and glacier mass change (R = 0.49 > 0.25 in Figure 8e) was not consistent with the results estimated by MLR analysis (p = 0.775 > 0.05, 95% confidence interval).

In order to further estimate the effect of morphometric factors on glacier mass changes at a more local scale, the whole ice field was divided into the grids with a size of 10 × 10 pixels. Note that the factor of glacier area was not included in the analysis due to the same size of all meshes. Figure 9 exhibits the correlation between glacier elevation changes and five morphometric factors over entire meshes of the PIF. Obviously, the factor of aspect was not a statistically significant control on glacier elevation changes. As shown in Figure 9a, there was not any apparent difference of glacier elevation changes in different aspects. The glacier mass changes might be slightly correlated with surface slope and curvature. Figure 9c,d illustrates that the grid with smaller slope and curvature might be more possible to experience a larger value of glacier thinning or thickening. Figure 9e shows that the factor of altitude should be an effective control of glacier elevation changes. In particular, the relationships between altitude and glacier elevation changes are different between the glacier tongue and upper regions. Specifically, when altitude increases, surface thinning over glacier tongue areas (altitudes of 5400–5700 m) decreases quickly, but surface thickening increases slowly around the ELA and its upper areas (altitudes of 5700–6200 m). It is noteworthy that glacier grids with higher ice flow rates possibly experienced a relatively stable glacier surface elevation, whereas obvious glacier thinning and thickening were usually found over quasi-stagnant grids (lower ice flow rate), as demonstrated in Figure 9b. It infers that surface velocities might affect glacier elevation changes at local scales.

Generally, the relationships between glacier elevation changes and morphometric factors over individual glaciers of the PIF were consistent with that at more local scale. However, obvious correlation between ice flow velocity and glacier mass change which was masked by averaging over individual glaciers was found for entire meshes on the ice field. Overall, according to the results of statistical analysis at different scales, we found that the observed spatially variable glacier mass changes might be related with the morphometric factors of altitude and velocity. Interestingly, the altitude-dependent glacier mass changes were basically supported by the significant correlations between altitude and main climatic parameters such as air temperature, relative humidity, air pressure, and precipitation (see Table 1 in Huintjes et al. [14]).

5. Discussion

5.1. Morphometric Effects on Glacier Mass Balances in the TP and Its Surroundings

In this study, our results revealed that the glacier mass balance of the Puruogangri ice field is primarily affected by mean altitude and ice flow velocity, and may be slightly affected by slope. Note that the debris cover is not included in our analysis because it does not cover glaciers of the PIF [29]. Generally, the altitude dependency of glacier mass balance over the PIF was supported by the stronger surface melting at lower elevations, which was found for several glaciers in the Himalaya and Karakoram [1,2]. However, Berthier et al. [19] reported that glaciers in the Spiti/Lahaul (western Himalaya) experienced thickening or limited thinning below 4400 m. Moreover, the velocity-dependent glacier elevation changes confirm the effect of velocity on glacier mass changes in Bhutan, Himalaya [23]. Additionally, although surface slope was slightly correlated with glacier mass change on the PIF, it was a major topographic control of glacier mass balance in eastern Himalaya [21].

Glacier mass balances of the PIF were not related with the factors of size and aspect, which is generally in agreement with the results reported by Berthier et al. [19] and Scherler et al. [24]. However, the glacier mass balances of some other glaciers in the adjacent regions were dependent on these factors. Specifically, the mass balances of glaciers in eastern Himalaya were generally determined by the factor of size [21]. An aspect-dependent glacier mass change was also observed in eastern Himalaya [23] and West Kunlun [20]. Consequently, for other glaciers in the TP and its surroundings, the effects of morphometric factors on glacier mass balances are not always consistent with our results. Furthermore, morphometric effects on glacier mass balance are far from homogeneous over different glaciated regions. Additionally, in one glaciated region, glacier mass changes are commonly affected by several morphometric factors. For example, both aspect and velocity are important controls of glacier mass balances in the Bhutan, Himalaya [23].

Overall, glaciers in the vicinity do not always experience the same mass change behaviors due to these local morphometric effects even in the same climatic setting. Consequently, a high uncertainty might exist in region-wide glacier mass change results by using the area-weighted extrapolation method from limited observations over several study sites in the Tibetan Plateau and its surroundings [4,23,44].

5.2. Climatic Influences on Glacier Mass Balance of the PIF

As shown in Figure 2, there was a slight mass gain in the interior areas, whereas obvious mass loss occurred in most glacier termini regions. To evaluate the effects of climate changes on the observed patterns of glacier mass balance over the PIF, we analyzed the precipitation and air temperature records between 1980 and 2012 from the Tuotuohe (34°13′ N, 92°26′ E), Bange (31°22′ N, 90°01′ E), and Anduo (32°21′ N, 91°06′ E) meteorological stations, located approximately 300 km from the study area. These meteorological data were downloaded from the China Meteorological Data Sharing Service System (http://data.cma.cn/).

Changes in glacier mass balance should be mainly related to spatial and temporal variations in summer temperature (June, July, and August) and annual precipitation in the Tibetan Plateau and its surroundings, because glacier melting usually dominates on the summer ablation season and snowfall occurs during all seasons of the year [2,45]. Figure 10a exhibits that summer temperature at Tuotuohe station is lower than two other stations due to higher latitude. Moreover, for temporal changes in summer temperature, a clearly positive trend was found over these three stations between 1980 and 2012, which was generally consistent with continuous warming trend since early 20th century [30]. Therefore, significant mass losses at the glacier terminus regions of the PIF were probably due to continuous increase in summer temperature.

Annual precipitation over the PIF (central TP) was much more stable than that in the Pamir, Karakoram, and Himalaya, because the climate of this study area is dominated by continental climatic conditions [2]. Furthermore, the ice field lies close to the boundary between the monsoon-dominated south and the westerly dominated north. Thus the decrease in precipitation caused by weakened Indian monsoon might be compensated by the strengthened westerlies. Figure 10b shows that the annual precipitations of these three stations were relatively stable during the period of 1980–2012. This result is basically agreement with the precipitation trend from 1979 to 2010 (Figure 4a in Yao et al. [2]). Particularly, it is noted that there was a slight increase of annual precipitation since 2000 (Figure 10b), which might contribute to the glacier mass gain over the interior areas of the ice field.

5.3. Surge-Type Glacier

Surge-type glacier is characterized by the life cycle of an active phase and a quiescent phase [46]. During the active phase, the glacier experiences a rapid transfer of mass from the upper “reservoir” area to the lower “receiving” area, and high ice flow rates that increase by a factor of 10 to 100, whereas ice flow rates decline substantially and the reservoir area builds up again during the quiescent phase [47]. According to the changes in glacier terminus location and surface velocity, some surge-type glaciers have been identified in the Karakoram [48,49], Pamir plateau [50], and Kunlun Shan [40,51]. Here, elevation changes probably induced by a surging event were clearly detected over the glacier 5Z213E0012 of the PIF, in the central Tibetan Plateau. Over the 12-year period (2000–2012), this glacier was thickened by up to 80 m at the terminus while thinned by approximately 20 m along its central flow line, as shown in Figure 3d. Moreover, Figure 11 demonstrates that glacier terminus advanced in the almost same time period (2000–2011) as that of the elevation change estimation, and continued in the following period of 2011–2016. These changes in glacier surface elevation and terminus locations imply that a surging event might be fully developed between 2000 and 2016. In addition, an obvious glacier retreat was observed over this glacier from 1976 to 2000 (Figure 11), which is a typical feature of the quiescent phase. Consequently, the surging event of glacier 5Z213E0012 might have started in 2000 and lasted for about 16 years.

Figure 12 illustrates the surface velocities in 2009 along the central flow line of glacier 5Z213E0012. In general, the magnitudes (0.001–0.045 m·day⁻¹) of our observed glacier velocities are similar with that (0.005–0.052 m·day⁻¹) in 1996 [14]. This finding is not consistent with that of most identified surge-type glaciers, which the velocity over active phase was 10–100 times more than that during quiescent phase [47], such as Karayaylak glacier [50] and central Yulinchuan glacier [51]. However, it is noteworthy that a significant change in ice flow pattern was found on glacier 5Z213E0012 between 1996 and 2009. As shown in Figure 12, a relatively higher ice flow rate was occurred on glacier upper portion (3000–6000 m from terminus) in 2009, whereas glacier tongue areas (0–3000 m from terminus) experienced a more rapid movement in 1996 (see Figure 3 in Huintjes et al. [14]). This change in glacier surface velocity might be caused by the surging event; however, further work is required to understand it, such as more detailed investigations of glacier dynamic mechanisms and bedrock characteristics.

6. Conclusions and Future Work

Clear spatial variability of mass balances for individual glaciers, ranging from −0.0351 to +0.0194 m·w.e.·year⁻¹, was observed on the Puruogangri ice field on the central TP from 2000 to 2012 based on a comparison of the TanDEM-X DEM (2012) with the SRTM-X DEM (2000). The method of linear regression analysis was employed to investigate the impacts of local morphometric factors on glacier mass balances. The results revealed a significant correlation (R = 0.84, p < 0.001) between the glacier mass balance and mean altitude, whereas glacier size, aspect, and ice flow velocity were not statistically significant controls. Furthermore, in order to further estimate the effect of morphometric factors on glacier mass changes at a more local scale, the whole ice field was divided into the grids with a size of 10 × 10 pixels. Importantly, apart from the factor of altitude, surface velocity is correlated with glacier elevation change at a local scale.

In addition, for the entire ice field, an evident pattern of mass loss was observed in most glacier tongue regions, whereas a slight mass gain was detected in the interior areas. This overall spatial behavior of glacier mass changes might be attributed to the regional climatic variations associated with the continuous increases in summer temperature between 1980 and 2012 and slight increase of annual precipitation since 2000. Future work will be conducted for comprehensively evaluating the impacts of climatic changes and morphometric factors on glacier mass balance, in order to provide greater insight into glacier changes and their effects on regional water resources and ecosystems.