Long-Term Deceleration of Glaciers in the Eastern Nyainqentanglha Mountains, Southeastern Tibetan Plateau, Revealed from Landsat Images

: Temperate glaciers are very sensitive to variations in temperature and precipitation, and thus represent a good indicator of climate change. By exploiting complete Landsat archives during periods of 1988–1990, 2000–2002 and 2014–2016, we derived three velocity maps of the temperate glaciers on the eastern Nyainq ê ntanglha Mountains in southeastern Tibetan to reveal the long-term changes of glacier surface velocity. Our results show that all the investigated glaciers experienced deceleration, with rates of deceleration varying from 4.15% to 29.8% per decade during the period from 1988–1990 to 2014–2016, showing heterogeneous deceleration patterns. A signiﬁcant rise in temperature and an insigniﬁcant decrease in precipitation was found from the meteorological data of the nearby meteorological station. The region-wide deceleration of glaciers was, thus, attributed to the negative mass balance induced mainly by the rise in temperature. The averaged rates of deceleration for periods from 1988–1990 to 2000–2002 and from 2000–2002 to 2014–2016 are 20.97% and 22.02% per decade, respectively, indicating a nearly even speed of deceleration in velocities during study periods. The nearly even speed of deceleration in velocities and the accelerating mass loss trend from periods before the 2000s to after the 2000s highlighted the complexity of the relationship between mass balance and glacier velocity. m 22 m / y from 2013 to 2014. These coincide with the velocities of 26 m / y, 78 m / y, and 41 m / y between 2000 and 2002 and 13 m / y, 42 m / y, 26 m / y from 2014 to 2016. The mean velocities of the Glacier were 40 m / y and 31 m / y for 2000–2002 and 2014–2016 in this study, respectively, whereas mean velocities were reported to be 52 m / y and 24 m / y for 1999–2003 and 2013–2014. The mean velocities of were 42 m m 2000–2002 in this study, respectively, whereas mean velocities were 32 m and m 1999–2003 are due


Introduction
Temperate glaciers, accounting for 18.6% of the total glacial number and 22.2% of the total glacier area in China, occur mainly in the east section and on the south slope of the Himalayas, the Hengduan Mountains, and the mid-east section of the Nyainqentanglha Mountains [1]. These glaciers are characterized by high accumulation and strong ablation, and are, thus, very sensitive to changes in temperature and precipitation, representing a good indicator of climate change [1][2][3]. In recent decades, the global mean annual air temperature experienced a remarkable increase, and most of the glaciers in the world are undergoing retreat, with particularly significant retreating trends in alpine glaciers [4,5]. In the Tibetan Plateau, the greatest shrinkage of a glacier took place in the southeastern Tibetan Plateau, where a decrease rate of 48.2 m yr −1 in length, and a shrinkage rate of 0.57% yr −1 in the area were observed during the 1970s-2000s [6].
Glacier velocity, which is a combination of ice deformation, bed deformation, and glacier sliding, is an important parameter to better study the dynamics of glaciers and their interplay with climate changes in the region [7][8][9]. In situ observations serve as one of the most accurate methods for measuring glacier velocity, but the remote areas where glaciers develop have prevented frequent visitation by

Data
We used Landsat satellite images to derive the multi-annual glacier velocity of the study glaciers. Images of path 135, row 39, and path 136, row 39 were needed to cover the study area completely. In order to derive reliable and spatially complete annual glacier velocity map for year T, we selected all pairs of the form (T − 1; T), (T; T + 1), and (T − 1; T + 1) to form a velocity stack so that all velocity measured was centered around year T, as described in [22]. The repeat cycle of Landsat satellites is 16 days, and 23 cycles represent 368 days. Therefore, we processed all the image pairs that were separated by 368 − 16, 368, and 368 + 16 days for one year and 736 − 16, 736, and 736 + 16 days for 2 years in the time spans 1988-1990, 2000-2002, and 2014-2016, respectively, to obtain three annual velocities of these periods.

Data
We used Landsat satellite images to derive the multi-annual glacier velocity of the study glaciers. Images of path 135, row 39, and path 136, row 39 were needed to cover the study area completely. In order to derive reliable and spatially complete annual glacier velocity map for year T, we selected all pairs of the form (T − 1; T), (T; T + 1), and (T − 1; T + 1) to form a velocity stack so that all velocity measured was centered around year T, as described in [22]. The repeat cycle of Landsat satellites is 16 days, and 23 cycles represent 368 days. Therefore, we processed all the image pairs that were separated by 368 − 16, 368, and 368 + 16 days for one year and 736 − 16, 736, and 736 + 16 days for 2 years in the time spans 1988-1990, 2000-2002, and 2014-2016, respectively, to obtain three annual velocities of these periods. From 1988 to 1990, we used satellite images of the Landsat 5 Thematic Mapper (TM), with 30 m resolution. From 2000 to 2002, we used images of the Landsat 5 TM in the combination of the multispectral images of Landsat 7 Enhanced Thematic Mapper Plus (ETM+), the resolutions of which are also 30 m. From 2014 to 2016, we used images of the Landsat 8 Operational Land Imager (OLI) panchromatic band (Band 8) because of the high resolution of 15 m and a large dynamic range. All the images were level 1T products, which are already terrain corrected using ground control points (GCPs) and digital elevation models (DEMs) and are freely available from the United States Geological Survey (USGS) EarthExplorer (http://earthexplorer.usgs.gov).

Methods
We followed the method proposed in [22] to derive the glacier velocity on the eastern Nyainqêntanglha Mountains for nearly three decades. First, some preprocessing was performed on images if needed. Then, velocity fields were extracted from each pair of images using the orientation correlation method. Finally, all the velocity fields in the same periods were merged using a spatiotemporal median filter, for each x and y component of the velocity, respectively, and the uncertainty was evaluated. In this section, we will describe the methods and the parameters used in this study briefly; readers are referred to [22] for more detailed information if necessary.

Preprocessing of Images
We derived glacier velocity using image pairs selected from three periods of 1988-1990, 2000-2002 and 2014-2016, respectively. As the images we used were level 1T products, topographical distortion was corrected. The displacement extracted from two images will not have been influenced by topography and represents the real horizontal motion. First of all, some images that were not exactly georeferenced were first co-registered to a reference image with subpixel accuracy. Images were then enhanced to improve the effectiveness of the correlation algorithm. For the images of Landsat 5/7, we applied principal component analysis (PCA) transformation on the images of bands 4 and 5. PCA transformation is a statistical procedure that converts a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components, by using an orthogonal transformation [34]. The first component of the PCA transformation on the images of bands 4 and 5 was proven to perform most robustly over mountain glaciers [22]. For Landsat 8, we applied the feature tracking method on the panchromatic band directly because of the high resolution and large dynamic range.

Extracting Velocity Fields from Each Image Pair
For each pair of images, the orientation correlation method operating in the frequency domain was used to produce velocity estimates between the first image called the reference image and the second image called the search image because of its efficiency and good performance over mountain glaciers [35]. First, the orientation images of the reference image and search image were generated using image gradient information. For an image I, its orientation image is a complex image in which the x direction gradient of I is the real part and the y direction gradient of I is the imaginary part.
The orientation image f is then normalized to take only the orientation into account: where g x = ∂I ∂x and g y = ∂I ∂y are x direction and y direction gradients of intensity I, respectively. Then, a reference window and a search window centered around pixel (i, j) were selected in the orientation image of the reference image and search image, respectively. These two windows were transformed into the frequency domain and matched using the orientation correlation method. The correlation surface OC for each pixel (i, j) was then given by: where IFFT indicates the inverse fast Fourier transform, F r (u, v) and F s (u, v) are the fast fourier transforms of the reference window and search window, respectively, and * means the complex conjugate. The maximum of the orientation correlation surface indicates the movement of the glacier, and the maximum is interpolated to obtain subpixel displacements. Before image matching, parameters such as reference window size, search window size, and spacing should be determined according to the image attributes and the characteristics of study glaciers. For 30 m resolution images, we set the reference window γ r to 16 × 16 pixels for Landsat 5/7 images and 32 × 32 for Landsat 8 images, covering an area of 480 × 480 m on the ground, which is the same as used in [22]. As displacements of glaciers in our study area can reach more than 400 m/y, we set the search window to vary from 48 to 128 pixels depending on the period and the pixel size. Spacing between two correlation patches was set to 4 pixels for Landsat 8 images and 2 for Landsat 5/7 images to generate an image of 60 m resolution.

Spatiotemporal Median Filtering and Uncertainty Assessment
After correlating all the selected image pairs, we first identified and removed mismatches or outliers by setting a threshold of signal-to-noise ratio (SNR), which is the ratio between the correlation maximum and the average value in the search window [36]. We chose a threshold of 5 to ensure that most of the suspicious matches could be eliminated while the more reliable estimates were reserved [22,37]. Then, all the velocity fields in the same periods were merged using a median filter both in a spatial and temporal neighborhood, for each x and y component of the velocity. Nine closest neighbors were included in the median filter, which covers an area 180 × 180 m on the ground. In order to obtain a statistically significant median, we removed the velocity vectors if the number of velocity estimates used to calculate the median was smaller than 5 [22]. This step also helped to remove spatially isolated pixels. As the median is not sensitive to isolated outliers, the spatiotemporal median filter could not only filter out isolated outliers that were not removed by the SNR threshold but also remove the incoherent errors caused by orthorectification errors, shadows, and clouds. Finally, each velocity map was inspected carefully, and some erroneous estimates were removed manually. The normalized median absolute deviation (NMAD) is a robust statistical description for the standard deviation, and it has been used in other applications [38,39]. During the spatiotemporal median filtering step, NMAD at each velocity pixel (i, j), symbolized as σ(i, j), can be formulated as follows [22]: where T is a set of velocity estimates, V(i, j, t) is the t th velocity elements of the set T, and V(i, j) is the median velocity obtained from the set T. The uncertainty of the spatiotemporal median filter result over the period considered, which is the 95% confidence interval of the final velocity distribution and symbolized as t 95 , is assumed to follow the law below [22]: where N is the number of velocity estimates which are used in the spatiotemporal median filter to compute the median velocity. σ is the NMAD of the N velocity estimates. k and α are parameters to be determined. A linear relationship was obtained after applying a logarithm to this equation, and the parameters to be determined became p 0 and p 1 : We computed the t 95 for pixels in the stable areas where the true velocity was known to be null and to fit the relationship between t 95 , σ, and N in Equation (5) using a least-square regression. By extrapolating the relationship to glacier areas, we were able to calculate the t 95 of each component for the final velocity.

Results
We processed all Landsat pairs available between 1988-1990, 2000-2002, and 2014-2016 in the study area, and three glacier velocity maps were generated by mosaicking two different paths, as shown in Figure 2. For better visualization, the velocity maps were subsetted by the Randolph Glacier Inventory version 5.0 (RGI v 5.0) (http://www.glims.org/RGI/rgi50_dl.html) and were superimposed on the hil-shade image generated from SRTM DEM. In general, the glaciers flow in almost the same patterns in three velocity maps, except that there is a slight decrease of glacier velocity when checked carefully. The maximum velocity occurs on Glacier 8, which is also the largest glacier in the area. There are some voids in the velocity maps, which occur mainly on the accumulation area of the glacier. Some of them also occur in the ablation area of some glaciers. Voids occurring on the accumulation area of the glaciers are caused, mainly, by the failure of image correlation due to the absence of features in these snow-covered areas. Voids appearing in the ablation area of some glaciers are mainly due to lack of features and the deformation of surface features due to strong surface velocity gradients or high ablation. Another source of voids comes from manually removing erroneous estimates. These erroneous estimates mostly appear miscellaneously and discontinuously. Thus, we can remove them manually by carefully inspecting the directions and values of the velocity measurements.
To quantitatively analyze changes in the glacier velocities, velocity profiles of the 26 investigated glaciers along the manually delineated centerlines for three periods were extracted ( Figure 3). Though the lower portions of most glaciers are covered by debris, the velocity data indicate that glaciers in the study area are still active. It can be observed that profiles along the centerlines for three periods are almost in the same pattern, which proves the reliability of the deriving results.
From the velocity profile data, the maximum velocity along the centerline of each glacier is larger than 20 m/y, and the maximum velocity among the 26 glaciers investigated occurred at Glacier 22, with a magnitude of more than 420 m/y. Generally, surface velocities are lowest at the glacier fronts but increase with elevation and affected by topography, tributaries, and other factors. For example, the almost linear increases of surface motion were observed for Glaciers 3, 4, 5, 6, and 7, until a maximum was reached at the end of the profiles (Figure 3c-g). However, glacier velocity may also be affected by the large slope (Glacier 2, as shown in Figure 3b) and the convergence of the glacier's main trunk and the large tributaries (Glacier 10, 13, 17, 20, 22, and 24, as shown in Figure 3j,m,q,t,v,x), in which cases large fluctuations can be observed from the velocity profiles of the glaciers. Though not in all cases, the velocities at the glacier termini of lake-terminating glaciers, such as for Glacier 8, 14, and 23, are usually larger than land terminating glaciers, the front of which usually seems to be stable with no obvious movement [12,40,41], as shown in Figure 3h,n,w. The velocities also depend on the different developmental stages of the proglacial lakes [42,43]. A higher surface velocity near the snout of a lake-terminating glacier is primarily driven by enhanced local ice mass loss due to increase in buoyancy at the calving front [44,45]; increasing basal lubrication due to the intrusion of lake water in the glacial sub-system [46]; and changes in the surface gradient of glacier [47]. Note that in Figure  A continuous and remarkable deceleration from 1988-1990 to 2014-2016 can be easily observed from the velocity data of the three periods, as shown in Figure 3. It is worth noting that glacier velocity changes between the periods 1988-1990, 2000-2002, and 2014-2016 cannot be evaluated, as the velocity fields derived here are composites of all available image pairs acquired during the period to increase the accuracy of the flow field. The statistical data for each glacier conducted from the profiles along the centerlines are summarized in Table 1. The mean velocity of each period in Table 1 was obtained by averaging all the velocity measurements along with the profiles, as done in [21]. As the profiles cover different portions of the glaciers; the mean velocities between glaciers are not comparable. Mean velocities can be used to reflect the velocity variations between different periods when velocity measurements can be calculated at the same portions. But mean velocities fail to represent the change of a glacier's velocity if velocity measurements are obtained at different portions. For example, in Figure 3c,h,k,x, deceleration can be seen for Glaciers 3, 8, 11, and 24, but the mean velocities at the three periods fail to reveal the deceleration (Table 1). To quantitatively analyze changes in the glacier velocities, velocity profiles of the 26 investigated glaciers along the manually delineated centerlines for three periods were extracted (Figure 3). Though the lower portions of most glaciers are covered by debris, the velocity data indicate that

Uncertainty Assessment
The uncertainty of the final velocity is decreased with an increasing number of observations and is also influenced by the NMAD of the velocities used to merge the median [22]. The relationships between the 95% confidence intervals of the final velocity in stable areas as a function of , i.e., the numbers of points used to compute the median and the NMAD for the three periods of two paths, are shown in Figure 4. In Figure 4d-f,j-l, there is a strong, linear relationship between ( / ) and ( ) for > 5 ( ( ) < 0.7).     Thus, in this study, we calculated the rate of deceleration by using the common velocity measurements that could be simultaneously obtained at two different periods to reduce this bias. The resulting percentage of common measurements and rate of deceleration between different periods are also given in Table 1

Uncertainty Assessment
The uncertainty of the final velocity is decreased with an increasing number of observations and is also influenced by the NMAD of the velocities used to merge the median [22]. The relationships between the 95% confidence intervals t 95 of the final velocity in stable areas as a function of N, i.e., the numbers of points used to compute the median and the NMAD σ for the three periods of two paths, are shown in Figure 4. In Figure 4d-f,j-l, there is a strong, linear relationship between log(t 95 /σ) and log(N) for N > 5 (log(N) < 0.7).
The parameters of the regression between log(t 95 /σ) and log(N) for log(N) < 0.7 (N < 5) are listed in Table 2. Using this relationship, we computed the uncertainty t 95 (i.e., 95% confidence interval) as a function of σ and N for every pixel. The median uncertainties in stable areas and on glaciers for [1988][1989][1990][2000][2001][2002], and 2014-2016 are given in Table 3, which are comparable with the uncertainties of 2 m/y in stable areas and 4.4 m/y on glaciers as described in [22]. The resulting uncertainties along the centerline profiles are presented as error bars, as shown in Figure 3. Generally, uncertainties of periods 1988-1990 and 2000-2002 are comparable, but uncertainties of periods between 2014-2016 are observably smaller than those of 1988-1990 and 2000-2002. This is because the image quality and geolocation of Landsat 8 are remarkably improved compared to Landsat 5/7 [48][49][50], contributing to the larger number of available observations and smaller dispersions. However, the uncertainties are small enough to characterize the glacier change between the study periods. The parameters of the regression between ( / ) and ( ) for ( ) < 0.7 ( < 5) are listed in Table 2. Using this relationship, we computed the uncertainty (i.e., 95% confidence interval) as a function of and for every pixel. The median uncertainties in stable areas and on glaciers for [1988][1989][1990][2000][2001][2002], and 2014-2016 are given in Table 3, which are comparable with the uncertainties of 2 m/y in stable areas and 4.4 m/y on glaciers as described in [22]. The resulting uncertainties along the centerline profiles are presented as error bars, as shown in Figure 3. Generally,

Comparison with Previous Studies
Studies of glacial velocity in the southeastern Tibetan Plateau are limited. Neckel et al. (2017) studied five glaciers that were covered by debris in the central-northern part of the Nyainqêntanglha Mountains from 1999-2014 [21]. These five glaciers correspond to Glaciers 14,15,17,18,19 in this study; thus, we compared our results of these five glaciers with the results in [21].  [21], as well as in this study. Considering the different times, center profiles, and different numbers of measurements obtained, our results proved to be consistent with those in [21].
In [21], Glacier 18 was found to experience the highest rate of deceleration of the five glaciers However, Neckel et al. (2017) calculated the rate of deceleration by using the mean velocities obtained from all the velocity measurements at two different periods [21]. In this study, we calculated the rate of deceleration by using the common velocity measurements that could be simultaneously obtained at two different periods to reduce this bias. By doing so, rates of deceleration of Glacier 15 and Glacier 18 between the period of 2000-2002 and the period of 2014-2016 were 48% and 54%, respectively. Among the five glaciers that were previously studied by Neckel et al. (2017), Glacier 18 still suffered the highest rate of deceleration between the period of 2000-2002 period and the period of 2014-2016, but Glacier 15 was not the one with the lowest deceleration rates. Instead, the lowest rate of deceleration was 23%, which corresponds to Glacier 14. The results in this study are also in general agreement with other glaciers' decelerations in many other places, such as the 43% per decade in Pamir, 25% per decade at the Penny Ice Cap, and 11% per decade in the Alaskan Range [51,52], and they in line with the neighboring Hailuogou Glacier, whose velocity decreased by 24% in ablation area from 1981 and 2008 [53].

Possible Reasons for Glacier Deceleration
Changes in mass balance redistributes ice mass to reach a new equilibrium influence, which would accordingly influence the ice velocities [54]. When the glacier is experiencing negative mass balance, its ice thickness and ice mass needed to transport down the glacier will reduce, resulting in reduced ice deformation. Ice velocities should reduce as a result of reduced ice deformation [7,55]. Glacier deceleration has also been reported in many places, such as the Pamir and Caucasus mountains, Penny Ice Cap, Alaskan Range, Patagonia, and southeast Greenland [52,56], and negative mass balance is considered to be the primary cause for the deceleration.
Over the past decades, air temperatures on the Tibetan Plateau have been seen to increase in many studies [57][58][59]. Due to the shortage of meteorological stations in the study area, data recorded at Bomi station (located at 29.528 N, 95.468 E, 2730 m a.s.l.) were collected to reveal the changing trend of temperature and precipitation during the study periods. The Bomi station is located 90 km west of our study area and is the nearest station. Since the Bomi station is also influenced by the South Asian monsoon, meteorological data at the Bomi station can be used to assess the trend of climate change in the study area. The mean annual temperature at the Bomi station increased significantly (p < 0.0001, R = 0.74, red dash line in Figure 5) during 1988-2016, at a rate of 0.40 • C per decade. While precipitation fluctuated, it witnessed a slight decrease, not statistically significant at the 5% level, from 1988-2016 (blue dash line in Figure 5).
respectively. Among the five glaciers that were previously studied by Neckel et al. (2017), Glacier 18 still suffered the highest rate of deceleration between the period of 2000-2002 period and the period of 2014-2016, but Glacier 15 was not the one with the lowest deceleration rates. Instead, the lowest rate of deceleration was 23%, which corresponds to Glacier 14. The results in this study are also in general agreement with other glaciers' decelerations in many other places, such as the 43% per decade in Pamir, 25% per decade at the Penny Ice Cap, and 11% per decade in the Alaskan Range [51,52], and they in line with the neighboring Hailuogou Glacier, whose velocity decreased by 24% in ablation area from 1981 and 2008 [53].

Possible Reasons for Glacier Deceleration
Changes in mass balance redistributes ice mass to reach a new equilibrium influence, which would accordingly influence the ice velocities [54]. When the glacier is experiencing negative mass balance, its ice thickness and ice mass needed to transport down the glacier will reduce, resulting in reduced ice deformation. Ice velocities should reduce as a result of reduced ice deformation [7,55]. Glacier deceleration has also been reported in many places, such as the Pamir and Caucasus mountains, Penny Ice Cap, Alaskan Range, Patagonia, and southeast Greenland [52,56], and negative mass balance is considered to be the primary cause for the deceleration.
Over the past decades, air temperatures on the Tibetan Plateau have been seen to increase in many studies [57][58][59]. Due to the shortage of meteorological stations in the study area, data recorded at Bomi station (located at 29.528 N, 95.468 E, 2730 m a.s.l.) were collected to reveal the changing trend of temperature and precipitation during the study periods. The Bomi station is located 90 km west of our study area and is the nearest station. Since the Bomi station is also influenced by the South Asian monsoon, meteorological data at the Bomi station can be used to assess the trend of climate change in the study area. The mean annual temperature at the Bomi station increased significantly (p < 0.0001, R = 0.74, red dash line in Figure 5) during 1988-2016, at a rate of 0.40 °C per decade. While precipitation fluctuated, it witnessed a slight decrease, not statistically significant at the 5% level, from 1988-2016 (blue dash line in Figure 5). The rise in temperature may decrease the accumulation area of the glacier, resulting in a reduction in snow accumulation. Therefore, ice flux, ice thickness, and hence, surface velocities, suffer a decrease [60]. Recent researches have shown that glaciers in the southeastern Tibetan Plateau are experiencing mass loss. By measuring the differences between the DEM generated from the Hexagon KH-9 and the SRTM DEMs, Zhou et al. (2018) showed that mass balance at the Yigong Zangbo basin from the mid-1970s to 2000 was −0.11 ± 0.14 m w.e.a −1 [61]. Neckel et al. employed the ICESat data to measure glacier mass changes on the Tibetan Plateau from 2003 to 2009, and they found the mass loss The rise in temperature may decrease the accumulation area of the glacier, resulting in a reduction in snow accumulation. Therefore, ice flux, ice thickness, and hence, surface velocities, suffer a decrease [60]. Recent researches have shown that glaciers in the southeastern Tibetan Plateau are experiencing mass loss. By measuring the differences between the DEM generated from the Hexagon KH-9 and the SRTM DEMs, Zhou et al. (2018) showed that mass balance at the Yigong Zangbo basin from the mid-1970s to 2000 was −0.11 ± 0.14 m w.e.a −1 [61]. Neckel et al. employed the ICESat data to measure glacier mass changes on the Tibetan Plateau from 2003 to 2009, and they found the mass loss in the eastern Nyainqentanglha range (−0.69 ± 0.36 m w.e.a −1 ) to be one of the highest [62]. Yao et al. (2012) showed that the mass balance of Parlung glaciers, which are in southeastern Tibetan Plateau, is nearly −1.1 m w.e.a −1 [6]. By fitting pixel-based linear elevation trends over ASTER DEMs acquired between 2000 and 2016, Burn et al. (2017) estimated the mass balance of Nyainqentanglha range to be −0.62 ± 0.23 m w.e.a −1 [63]. Though not exactly at the same time, these previous studies have highlighted the fact that the Nyainqentanglha range is losing mass, with an accelerating mass loss trend from periods before the 2000s to periods after the 2000s. Thus, the region-wide deceleration of glaciers in the study area should be attributed to the negative mass balance that is mainly caused by the rise in temperature in this area. However, the averaged rates of deceleration for periods from [1988][1989][1990]

Conclusions
Temperate glaciers, which account for nearly one-quarter of the total glacial area in China, and are mainly in the southeastern Tibetan Plateau, are very sensitive to variations in temperature and precipitation. Under the influence of global warming, glaciers in the southeastern Tibetan Plateau have shrunk heavily. Observations of glacier velocities can be helpful for a better understanding of the dynamics of glaciers and their relationships to climate changes.
In this paper, to reveal the long-term changes of glacier surface velocity in the eastern Nyainqêntanglha Mountains, we derived velocity maps of [1988][1989][1990][2000][2001][2002], and 2014-2016 by exploiting complete Landsat archives during those periods, and the uncertainty was analyzed. The results show that glaciers are still in active state in the eastern Nyainqêntanglha Mountains. All the glaciers that have been studied experienced deceleration but with heterogeneous rates of deceleration, varying from 10.78% to 77.46% during the period from 1988-1990 to 2014-2016, or 4.15% and 29.8% per decade. The results are in general agreement with previous results in this area and other regions of High Mountain Asia. A significant rise in temperature and an insignificant decrease in precipitation can be found from the meteorological data of the nearby meteorological station and previous studies. The deceleration was, thus, considered to be induced by the negative mass balance, which mainly resulted from the rise in temperature. A nearly even speed of deceleration in glacier velocities was observed, as the averaged rates of deceleration for periods from 1988-1990 to 2000-2002 and from 2000-2002 to 2014-2016 were 25.17% and 30.83%, or 20.97% and 22.02% per decade, respectively. The nearly even speed of deceleration in glacier velocities together with an accelerating mass loss trend from periods before the 2000s to after the 2000s revealed by previous studies, highlighted the complex relationship between glacier velocity and mass balance, which can also be influenced by the local topographic, climatic context and subglacial hydrology, amongst other factors [64,65]. The heterogeneous deceleration patterns of individual glaciers in this area are also a combined result of these factors functioning together.
Author Contributions: X.F. designed the research, performed the experiments, and wrote the paper. J.Z. supervised the work and revised the manuscript.