Multitemporal Glacier Mass Balance and Area Changes in the Puruogangri Ice Field during 1975–2021 Based on Multisource Satellite Observations

: Due to climate warming, the glaciers of the Tibetan Plateau have experienced rapid mass loss and the patterns of glacier changes have exhibited high spatiotemporal heterogeneity, especially in interior areas. As the largest ice ﬁeld within the Tibetan Plateau, the Puruogangri Ice Field has attracted


Introduction
The cryosphere is one of the most sensitive geospheres to global climate change and glaciers especially, as important features of the cryosphere, have always been regarded as an indicator of climate change [1][2][3].Global glaciers (including ice sheets) cover an area of 14.51 × 10 6 km 2 and the volume of total ice reserves is about 27.6 × 10 6 km 3 , which stores 75% of global freshwater resources [4].The Tibetan Plateau (TP) is the largest modern glacier area in the middle and lower latitudes of the world and the glaciers on most parts of the TP have generally been undergoing an accelerated retreat due to climate warming [3][4][5], which has had a great impact on water resources [6], sea-level rise [7] and increased glacial hazards, such as glacier collapse [8,9] and glacial lake outburst floods [10,11].By studying the changes in mountain glaciers on the TP, it is not only possible to understand the characteristics and trends of glacier changes in this region but also to further investigate the relationships between glacier changes and climate change, regional water availability and disaster occurrence, which is of great importance [12][13][14].Recent studies on glacier changes have mainly focused on changes in glacier length, area, flow velocity, thickness and mass balance [12,15].The few earlier accessible glacier studies were mainly based on field measurements and their results were generally regarded as more reliable [16].However, as glaciers are usually located in remote areas on the TP, it has been difficult to conduct large-scale and long-term in situ observations [3,17] due to the high costs of human, material and financial resources.With the development of remote sensing (RS) technology, satellite images have been increasingly widely used for the study of glacier changes and their measurement accuracy and efficiency have been continuously improved [18,19].Recent studies have shown that the TP glaciers, overall, have been in a state of pronounced mass loss and that the patterns of glacier changes have been highly heterogeneous at spatiotemporal scales [20].The rates of glacier changes have been inconsistent, especially in the interior of the TP (which is characterized by a relatively sparse distribution of glaciers), possibly due to different climatic responses in local areas [2,14,21].Hence, it is of great significance to carefully monitor glacier changes in detail over long periods and multiple time scales at main glacier centers to improve our understanding of the responses and adaptations of glaciers to local climate change.
The Puruogangri Ice Field and its surroundings, as the largest glacier center within the TP, is an ideal natural testing ground for studying glacial changes.For example, Zhou et al. [22] used the band ratio and manual interpretation methods to extract data for glacier area parameters in the Geladaindong mountain region and found that the glacier area decreased by 66.68 km 2 (7.37%) from 1992 to 2009 and shrank by 19.90 km 2 (2.84 km 2 a −1 ), 39.95 km 2 (5.71 km 2 a −1 ) and 6.83 km 2 (2.28 km 2 a −1 ) over the three time periods of 1992-1999, 1999-2006 and 2006-2009, respectively.Ye et al. [23] reported that glaciers in the Geladaindong mountain region had shrunk at a rate of 1.29 km 2 a −1 during 1969-2002, which was derived from Landsat images and topographic maps.As an important research parameter, glacier mass balance is a direct metric of glacier change and an important factor for exploring the response of glaciers to climate change.Research methods for glacier mass balance mainly include model estimation [13] and RS monitoring [13,14,24], especially the widely applied DEM differential method.Zhang et al. [25] derived the annual mass balance of three glaciers in the interior of the TP using an albedo-based method and found that the mass loss rate over 2000-2016 was 0.535 ± 0.063 m w.e. a −1 , 0.243 ± 0.066 m w.e. a −1 and 0.113 ± 0.068 m w.e. a −1 for the Xiao Dongkemadi glacier, Geladaindong mountain region glaciers and Puruogangri Ice Field, respectively.Studies that are based on RS data, such as satellite altimetry data and multisource DEMs, have been used to estimate mass balance by obtaining information about glacier surface elevation changes.Zhou et al. [26] obtained a glacier mass balance result of −0.22 ± 0.12 m w.e. a −1 for the Geladaindong region during 1970s-2000s, based on KH-9 images and SRTM DEM.Liu et al. [27] obtained a mean elevation change of about −0.14 ± 0.26 m a −1 for the Geladaindong mountain glaciers over 2000-2014 using SAR interferometry on TanDEM-X data.
The findings of the above area and mass balance studies have suggested that glaciers in the central TP have suffered from continuous mass loss.As the largest ice field within the central part of the TP [28], measuring the area change and mass balance of the Puruogangri Ice Field over long time scales could effectively reflect the overall trends of glacier changes in the central TP that are caused by climate change.Neckel et al. [29] reported slight mass losses of −0.038 ± 0.023 m w.e. a −1 and −0.044 ± 0.015 m w.e. a −1 in the Puruogangri Ice Field during 2000-2012 using the common DEM differencing method and the differential interferometric synthetic aperture radar (DInSAR) technique.Liu et al. [30] estimated an annual surface thinning rate of −0.317 ± 0.027 m a −1 between 2012 and 2016 by differencing TanDEM-X DEMs.By differencing five pairs of TanDEM-X DEMs that were acquired between April 2012 and January 2016, Liu et al. [21] reported a mass gain of 0.44 ± 0.10 m w.e. a −1 over 2011-2012 and mass losses of −0.13 ± 0.03, −0.34 ± 0.06 and −0.52 ± 0.10 m w.e. a −1 over 2012-2013, 2013-2014 and 2014-2016, respectively.However, the above-mentioned works mainly focused on the period of 2000-2016 and there is still a need to use RS data and methods for observations from before 2000, as well as more recent years.Therefore, conducting studies on pre-2000 periods and more recent periods could provide a more precise understanding of specific changes in the Puruogangri Ice Field and help to understand its responses to long-term climate change.
The main objective of this study was to update the area changes in the Puruogangri Ice Field and estimate the glacier mass balance over different periods from 1975 to 2021 using a higher temporal resolution than that in existing studies.The first step was to calculate the area changes for 1975-2021 using the band ratio method with Landsat and KH-9 images.Secondly, we quantified decade-scale changes in glacier thickness and mass balance for 1975-2000 and 2000-2021 using geodetic methods with KH-9 DEM, SRTM DEM and ZY-3 DEMs.On the basis of the previous point, we carried out intensive observations of the mass balance over 2000-2021 using Copernicus DEM and ZY-3 DEMs.Thirdly, observations from nearby meteorological stations were analyzed and their effects on glacier mass changes were examined.

Study Area
The Puruogangri Ice Field (33 • 43 ~34 • 04 N, 89 • 00 ~89 • 20 E) lies in the central part of the TP to the northeast of Shuanghu County (Figure 1) and has an area of 398.25 km 2 , according to the Randolph Glacier Inventory, version 6.0 (RGI 6.0) [29].It is the largest ice field within the TP and is on one of the highest mountains in the central area, with an average elevation of 5820 m a.s.l.(with the highest point at 6482 m a.s.l. and the lowest point at 5350 m a.s.l.).The Puruogangri Ice Field is generally flat with a northwestsoutheast strip distribution and emits more than 50 glacier tongues of varying lengths into the surrounding valleys.
The climate in this region is dominated by a cold and dry continental climate [28,31].According to the records from meteorological stations, precipitation in the region is scarce, with an annual precipitation level of about 100~500 mm.The high altitude makes the temperature generally low, with an average temperature of less than 1.0 • C. The records from two ice cores that were collected in 2000 and meteorological stations have indicated a significant warming since the 20th century and the shrinking range of the Puruogangri Ice Field has also been increasing [21,31].The climate in this region is dominated by a cold and dry continental climate [28,31].According to the records from meteorological stations, precipitation in the region is scarce, with an annual precipitation level of about 100~500 mm.The high altitude makes the temperature generally low, with an average temperature of less than 1.0 °C.The records from two ice cores that were collected in 2000 and meteorological stations have indicated a significant warming since the 20th century and the shrinking range of the Puruogangri Ice Field has also been increasing [21,31].

Data
The data that were adopted in this study are summarized in Table 1, in which brief descriptions of these data are presented.

Data
The data that were adopted in this study are summarized in Table 1, in which brief descriptions of these data are presented.The Landsat series (which are jointly managed by NASA and the USGS) have launched eight satellites so far since 1972 that have been configured with multispectral scanner (MS), thematic mapper (TM), enhanced thematic mapper plus (ETM+) and operational land imager (OLI) sensors, which provide continuous and effective data support for the longterm monitoring of glacier changes.The quality of Landsat images and the spatial resolution of their optical bands have been gradually improved in recent satellites, e.g., the highest resolution has been increased to 15 m in the Landsat 7/8 panchromatic band.In this study, we utilized one Landsat TM image from 2000, one Landsat ETM+ image from 2012 and three Landsat OLI images from 2015, 2017 and 2021 to extract the glacier outlines within the Puruogangri Ice Field so we could monitor the area changes (Table 1).The Landsat L1TP (Level 1 Precision Terrain) images were obtained from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/,accessed on 26 October 2021).
(2) KH-9 The KH-9 Hexagon is a photographic reconnaissance satellite program that includes a total of 20 series satellites, which were launched between 1971 and 1986 and have a high resolution of 6-9 m and a single scene coverage of approximately of 3400 km 2 that enable large area stereoscopic observations [32].In 2002, the USGS declassified about 29,000 images that were taken by the KH-9 satellite mission [33].In this study, we delineated the outlines and extracted DEMs for the Puruogangri Ice Field using a pair of KH-9 stereo images from 20 December 1975 (Table 1), which were downloaded from the USGS (https://earthexplorer.usgs.gov/,accessed on 2 October 2021).
(3) ZY-3 The ZY-3 is China's first civilian multispectral satellite, which was launched on 9 January 2012 and is equipped with three panchromatic cameras and one multispectral camera.The resolution of the panchromatic cameras is 3.5 m in forward (FWD) and backward (BWD) visuals and 2.1 m in nadir (NAD) visual, which meet the requirements for higher stereo mapping precision [34].The base to height ratio (B/H) of ZY-3 stereo images is 0.89, which is sensible for DEM generation [35].In this study, we purchased three ZY-3 stereo-pair images to extract DEMs for the years of 2015, 2017 and 2021 (Table 1).

DEM Products (1) SRTM DEMs
The Shuttle Radar Topography Mission (SRTM) was completed on 22 February 2000 after an 11-day mission and provided a high-quality DEM dataset for almost all global land area between 60 • N and 57 • S for the first time [36,37].In this study, the 1 arc-second SRTM C-band DEM was utilized to acquire the glacier surface elevation in 2000 and the X-band DEM was used to estimate the penetration depth of the C-band [26] (Table 1).It is worth noting that the datum references for the C-and X-band DEMs were different.The horizontal and vertical references of the C-band DEM were the WGS84 ellipsoid and EGM (Earth Gravitational Model) 1996 geoid, respectively, while those of the X-band DEM were both WGS84 ellipsoid [12].The SRTM C-and X-band DEMs were acquired from the USGS (https://earthexplorer.usgs.gov/,accessed on 19 August 2022).
(2) Copernicus DEM The Copernicus DEM is a digital surface model (DSM) that represents the Earth's surface, including buildings, infrastructure and vegetation [38].The Copernicus DEM is a WorldDEM product that was generated using radar satellite data that were acquired from the TanDEM-X mission [38,39], which provided the most realistic representation of height and should be one of the most important choices for quantitative topographic analyses [40,41].The Copernicus DEM provides three different products, which are named EEA-10, GLO-30 and GLO-90, and we used Copernicus DEM GLO-30 in this study.The original TanDEM-X images that covered the study area were acquired around July 2012.We did not account for the effects of X-band penetration in July that was due to summer surface melting [42,43].The horizontal and vertical datum references for the Copernicus DEM were WGS84 and EGM2008, respectively [38,44].The data were downloaded from the FTP server (https://panda.copernicus.eu/,accessed on 29 November 2021).

Satellite Altimetry Data
The ICESat-1 (Ice, Cloud and Land Elevation Satellite) was launched by NASA on 13 January 2013 and was the world's first Earth observation laser altimetry satellite with a revisit period of 91 days [45,46].The main payload is GLAS (the Geoscience Laser Altimeter System) and the footprint and spacing are about 70 m in diameter and roughly 172 m along the track direction, respectively [18,47,48].In this study, we used the GLAH14 product (global land surface altimetry data) from 2003 to 2009, which has increased Gaussian elevation control points to improve its accuracy, to assess the accuracy of the DEMs [18].The GLAH14 product was downloaded from the National Snow and Ice Data Center (NSIDC) (https://nsidc.org/data/GLAH14/versions/34,accessed on 23 October 2021).

Methods
The workflow that was followed to monitor the glacier area and mass balance changes in the Puruogangri Ice Field using multisource RS data is illustrated in Figure 2. As well as the SRTM and Copernicus DEMs, the DEMs for 1975 and three recent periods (i.e., 2015, 2017 and 2021) were extracted from KH-9 and ZY-3 stereo image pairs, respectively (Table 1).The above-mentioned DEMs were then used to estimate the glacier elevation changes over different periods by differentiating the DEMs of different years and eliminating the effects of planimetric and altimetric biases.Landsat images for 2000, 2012, 2015, 2017 and 2021 were used to delineate the corresponding boundaries using the band ratio method and the glacier boundaries from 1975 were delineated from a KH-9 image using visual interpretation.These were then used to estimate the changes in glacier area over different periods.Finally, the changes in glacier mass balance were estimated using the glacier area, elevation and density parameters from the different periods.Detailed descriptions of the key methods and processes are presented below.

KH-9 DEM
A common method for extracting DEMs from a pair of KH-9 stereo images is photogrammetry, which generally includes three steps.In this study, distortion correction and image splicing were performed first, then ground control points (GCPs) were manually selected in the ERDAS (Earth Resource Data Analysis System) LPS (Leica Photogrammetry Suite) and tie points (TPs) were automatically generated.Then, epipolar images were obtained and finally, the DEM was obtained after post-processing [14].It should be noted that the horizontal and vertical references that were used for the collection of the GCPs were derived from Google Earth images from various years and large differences in the spatiotemporal resolution of the satellite images increased the difficulty in collecting the GCPs.In our study, we adopted the automatic methodology that was presented by Maurer et al. [32], which uses a combination of traditional photogrammetry and computer vision with no requirement for the manual collection of GCPs.In the first instance, image stitching and cross positioning were completed using automatic feature matching and geometric transformation.Local adaptive filters, noise filters and histogram matching were applied to enhance the images after the separation of the processing grid windows.Secondly, the SURF (speeded up robust features) algorithm was used to extract the feature points from the stereo images and then those feature points were automatically matched.The SGBM (semiglobal block-matching) algorithm was used for the stereo matching of the left-and right-hand images.Finally, the point cloud coordinates were converted into the world coordinate system and co-registered to the SRTM C-band DEM.Due to the low optical contrast at high altitudes, the generated DEM generally showed high accuracy over flat areas but presented voids on steep slopes.

KH-9 DEM
A common method for extracting DEMs from a pair of KH-9 stereo images is photogrammetry, which generally includes three steps.In this study, distortion correction and image splicing were performed first, then ground control points (GCPs) were manually selected in the ERDAS (Earth Resource Data Analysis System) LPS (Leica Photogrammetry Suite) and tie points (TPs) were automatically generated.Then, epipolar images were obtained and finally, the DEM was obtained after post-processing [14].It should be noted that the horizontal and vertical references that were used for the collection of the GCPs were derived from Google Earth images from various years and large differences in the spatiotemporal resolution of the satellite images increased the difficulty in collecting the GCPs.In our study, we adopted the automatic methodology that was presented by Maurer et al. [32], which uses a combination of traditional photogrammetry and computer vision with no requirement for the manual collection of GCPs.In the first instance, image stitching and cross positioning were completed using automatic feature matching and geometric transformation.Local adaptive filters, noise filters and histogram matching were applied to enhance the images after the separation of the processing grid windows.Secondly, the SURF (speeded up robust features) algorithm was used to extract the feature

Glacier Boundary Delineation
The glacier boundaries within the Puruogangri Ice Field were extracted using the semi-automated band ratio method (Red/SWIR) [18] with band 3 and band 5 of the Landsat TM and ETM+ images from 2000 and 2012 and with band 4 and band 6 of the Landsat OLI images from 2015, 2017 and 2021.The first step was radiation correction, which included radiation calibration and atmospheric correction [10].Then, the band ratio method was employed using a threshold of 2.0 to distinguish glaciers from non-glacial terrain.The obtained boundaries were jagged as a result of the properties of the Landsat images and needed to be further smoothed.Finally, further visual interpretation was required to match the actual glaciers using Landsat false-color images.The glacier boundaries for 1975 were delineated from a KH-9 image using visual interpretation and were adjusted manually according to the glacier boundaries in 2000.

C-Band and X-Band Penetration Depth
The penetration depth of a SRTM C-band radar beam in ice and snow could result in lower elevation values over glacier surfaces, which had to be considered [10,12,52].Since the X-band penetration was much smaller than that of the C-band, we estimated the penetration depth of the C-band radar beam by differentiating the SRTM X-band DEM and C-band DEM on glacier surfaces [19].The first step was to convert the geographic coordinates of the SRTM C-band and X-band DEMs into projection coordinates and resample them to a resolution of 30 m.In addition, the X-band DEM needed to be converted into EGM96 to unify the height datum [53].Then, the biases of the different images in the X-and C-band DEMs were corrected using the processing method with elevation difference maps.
The C-band penetration depth was calculated for each 100-m altitude band of the Puruogangri Ice Field (Figure 3).The penetration depth in the ice-free areas was close to zero, while it ranged from 0.40 m to 5.39 m in the glaciers.The average penetration depth for the ice-free and glacier regions were 0.26 ± 0.02 and 2.25 ± 0.13 m, respectively, which were in agreement with the result of 2.2 m that was found by Liu et al. [30].Zhou et al. [26] measured the penetration depth in the Geladaindong region as 1.9 m by differencing SRTM-C and SRTM-X DEMs and Li et al. [54] obtained a value of 2.25 ± 0.09 m in the central TP using the same method.These findings also suggested that our results were reasonable.Although the X-band radar penetration depth was also a possible source of systematic bias for glacier mass balance, we did not consider the differences in the penetration depth of the X-band radar in glacier surfaces for the Copernicus DEM (which was acquired in July) because of the presence of meltwater in the summer.For the unconsidered X-band penetration depth in previous penetration estimates for the SRTM C-band DEM, we adopted the X-band penetration depth that was reported by Liu et al. [21], which was 0.61 ± 0.06 m in the study region.The DEM differencing method was employed in this study to monitor glacier elevation changes in the Puruogangri Ice Field using various DEMs from different periods be-

Estimation of Glacier Elevation Changes
The DEM differencing method was employed in this study to monitor glacier elevation changes in the Puruogangri Ice Field using various DEMs from different periods between 1975 and 2021, i.e., a KH-9 DEM from 1975, an SRTM-C DEM from 2000, a Copernicus DEM from 2012 and DEMs that were extracted from ZY-3 images from 2015, 2017 and 2021.The three-step method that was proposed by Nuth and Kääb [55] was used to co-register and improve the deviations between the DEMs from different years, which was based on differences in elevation, slope and aspect over stable terrain.The reason for adopting this three-step method was that there were three potential biases, including bias that was related to geographic positioning, elevation-dependent bias and deviations that were associated with data acquisition geometry [55].The elevation bias was divided into horizontal and vertical shift vectors and the employed analytic equation was: where ∆h 1 is the elevation difference over stable terrain, a is the horizontal shift vector, b is the direction of the horizontal shift, β is the aspect, α is the slope and (dh) is the vertical shift of the system [10].The equation was solved using the least-squares method, according to the elevation difference maps of the stable terrain.Then, the elevation-dependent bias was removed using a linear or polynomial relationship.Figure 4 shows the distribution of the normalized elevation differences between the SRTM-C DEM from 2000 and the Copernicus DEM from 2012 before and after co-registration, as an example.It can be clearly observed that after co-registration, the functional relationship between the differences in elevation, slope and aspect significantly improved.
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 After co-registration, planimetric trend bias and altimetric bias remained, w were associated with the position and accuracy of the GCPs during DEM genera [26,56].In this study, the planimetric trend bias (∆ℎ ) was corrected using quadratic face fitting (Equation ( 2)) and the altimetric bias (∆ℎ ) was corrected using quadratic ynomial simulation (Equation (3)): where  and  denote the horizontal and vertical coordinates of an image pixel p respectively,  is the value of an image pixel, ∆  and ∆  are the elevation differe at that point and   and   represent the model parameters to be estimated (the le squares method was used to solve these parameters).Figure 5 shows the histogram tribution of elevation differences over stable regions during 1975-2000 before and planimetric and altimetric corrections, as an example.After co-registration, planimetric trend bias and altimetric bias remained, which were associated with the position and accuracy of the GCPs during DEM generation [26,56].In this study, the planimetric trend bias (∆h 2 ) was corrected using quadratic surface fitting (Equation ( 2)) and the altimetric bias (∆h 3 ) was corrected using quadratic polynomial simulation (Equation (3)): (2) where x and y denote the horizontal and vertical coordinates of an image pixel point, respectively, z is the value of an image pixel, ∆h 2 and ∆h 3 are the elevation differences at that point and a i and b i represent the model parameters to be estimated (the least-squares method was used to solve these parameters).Figure 5 shows the histogram distribution of elevation differences over stable regions during 1975-2000 before and after planimetric and altimetric corrections, as an example.
where  and  denote the horizontal and vertical coordinates of an image pixel point, respectively,  is the value of an image pixel, ∆  and ∆  are the elevation differences at that point and   and   represent the model parameters to be estimated (the leastsquares method was used to solve these parameters).Figure 5 shows the histogram distribution of elevation differences over stable regions during 1975-2000 before and after planimetric and altimetric corrections, as an example.

Calculation of Glacier Mass Balance
The glacier volume changes were computed based on the changes in glacier area and thickness by multiplying the average glacier elevation change within each elevation interval (100 m) by the corresponding area and then summing the results [10,19,53].It is worth

Calculation of Glacier Mass Balance
The glacier volume changes were computed based on the changes in glacier area and thickness by multiplying the average glacier elevation change within each elevation interval (100 m) by the corresponding area and then summing the results [10,19,53].It is worth mentioning that it was necessary to fill in data voids due to data accuracy and the aforementioned data processing.Assuming that glacier elevation changes were similar over specified altitude ranges, we filled the data voids using the average elevation changes at those altitudes, i.e., the glacier was divided into altitude bands with 100-m intervals and the average elevation changes within each altitude band were calculated to fill the data voids [19,26].Then, a density of 850 ± 60 kg m −3 was adopted to convert the volume changes into glacier mass changes [57].The glacier mass balance was estimated using: where (∆h i ) is the average glacier elevation change for each altitude interval, A i is the glacier area for each altitude interval, n represents the number of altitude intervals, A represents the whole glacier area, ρ g is the glacier mass density and ρ w is the water density (1000 kg m −3 ).

Uncertainty Assessment
The uncertainty of the glacier boundary delineation was estimated using the length of the glacier boundaries and the image resolution [10]: where l is the length of a glacier boundary and r is the resolution of the satellite images.
For the uncertainty of the glacier elevation differences (δ ∆h ), we used the standard deviation of the elevation differences (δ E ) over stable regions and the autocorrelation distances (D).Then, the value of δ ∆h could be calculated using three formulae: ) where N e f f is the number of effective independent measurements, N t is the total number of measurements, R is the pixel resolution and D is the spatial autocorrelation distance.For the value of D, 600 m was adopted [56].Thus, the uncertainty of the mass balance (δ m0 ) was computed using the error propagation law (Equation ( 8)).In addition, since the uncertainty of the differences in the penetration depth between the C-band and X-band (δ C−X ) was 0.13 m and the uncertainty of the X-band penetration depth (δ X ) was 0.06 m, the final uncertainty of the mass balance (δ m ) was computed using Equation (9): ) where S is the glacier area, ∆h is the difference in glacier elevation, δ s and δ ρ g are the uncertainty of the area and glacier mass density, respectively, and S t is the total area.

Changes in Glacier Area
Table 2 shows the estimated glacier areas in the Puruogangri Ice Field for six different years.It was found that the total area of the Puruogangri Ice Field shrank from 1975 to 2017 and became stable between 2017 and 2021.On the whole, the glacier area shrank by −39.57± 1.41 km 2 during 1975-2021.Figure 6a shows the extracted boundaries for each year and Figure 6b shows the overall area changes for each period.From Figure 6a, we could see that from 1975 to 2000, the Puruogangri Ice Field generally showed continuous retreat, especially the G98 glacier in the north, the G94 glacier in the southwest and the G37 glacier in the southeast.After 2000, most of the glaciers continually shrank, while some glaciers advanced, such as the G05 glacier and the G13 glacier in the east.It could be observed that glaciers with significant area changes tended to have larger termini or were long and narrow in shape.As illustrated in Figure 6b, the overall area change rate during 1975-2021 was −0.86 km 2 a −1 , with rates of −0.83 km 2 a −1 and −0.90 km 2 a −1 for 1975-2000 and 2000-2021, respectively, which showed a slight increase in the shrinkage rate over the two most recent decades.More specifically, the rate of area change was −0.59 km 2 a −1 for 2000-2012 and −1.31 km 2 a −1 for 2012-2021 on a decade scale, which indicated a larger shrinkage over the past decade.For the period of 2012-2021, the rate at a multiyear scale was −2.38 km 2 a −1 , −2.74 km 2 a −1 and 0.40 km 2 a −1 for 2012-2015, 2015-2017 and 2017-2021, respectively.The rate of area change was the highest during 2015-2017, while the area expanded during 2017-2021 due to the surge of the G98 glacier (Figure 6).To better understand the spatial variations in the glacier area changes in the Puruogangri Ice Field, we selected six typical glacier termini (Figure 6a) for further investigation.The estimated areas of the selected glacier termini in six different years are also presented in Table 2 and the rates of area change in the different periods are shown in Figure 7.In general, the selected glacier termini showed three patterns of area changes: continuous retreat (i.e., the G94, G37 and G36 glaciers, as shown in Figure 7a-c); discontinuous retreat (i.e., the G05 glacier, as shown in 7d); and retreat and then advance (i.e., the G13 and G98 glaciers, as shown in Figure 7e, f)).As examples of the first pattern, the G94, G37 and G36 glaciers (Figure 7a-c) shrank continuously and decreased by −1.26 km 2 (−0.05 km 2 a −1 ), −3.83 km 2 (−0.15 km 2 a −1 ) and −0.73 km 2 (−0.03 km 2 a −1 ) between 1975 and 2000, respectively, and −2.20 km 2 (−0.10 km 2 a −1 ), −3.04 km 2 (−0.14 km 2 a −1 ) and −1.07 km 2 (−0.05 km 2 a −1 ) between 2000 and 2021, respectively.As an example of the second pattern, the area change of the G05 glacier (Figure 7d) was erratic and shrank −0.19 km 2 during 1975-2000, expanded 0.53 km 2 during 2000-2015 and then started to shrink again, decreasing from 14.62 ± 0.31 km 2 (2015) to 14.40 ± 0.31 km 2 (2021).As examples of the third pattern, the areas of both the G13 and G98 glaciers (Figure 7e, f) shrank from 1975 to 2017 with changing rates of −0.02 km 2 a −1 and −0.14 km 2 a −1 , respectively.Then, from 2017 to 2021, a surging phenomenon occurred in the G98 glacier, which led to an area increase of +1.34 km 2 (+0.34 km 2 a −1 ).The area of the G13 Glacier also slightly increased by +0.32 km 2 (+0.08 km 2 a −1 ) during this period.

Changes in Glacier Elevation and Mass Balance
Figure 8 shows the spatial distribution of the estimated glacier elevation changes in the Puruogangri Ice Field over different periods.The glaciers of the Puruogangri Ice Field thinned between 1975 and 2000, with significant surface reductions concentrated at the glacier termini, especially those of the G98, G36 and G37 glaciers (Figure 8a).Between 2000 and 2021 (Figure 8b), several glacier termini continued to thin, especially those of the G94 glacier, and showed more thinning over this period in comparison to that over 1975-2000.On the contrary, some thickening phenomena was observed at the termini of several other glaciers (e.g., the G98 and G05 glaciers).To describe the changes in the glaciers more clearly, Figure 9 shows the results of the overall glacier mass balance over different periods.It was found that the mass balance of the Puruogangri Ice Field was −0.23 ± 0.02 m w.e. a −1 and −0.29 ± 0.02 m w.e. a −1 during 1975-2000 and 2000-2021, respectively, which indicated a continuous and steady state of mass loss over the last four decades.
To further define the spatial heterogeneity of the glacier changes over the past 20 years, the elevation changes during the periods of 2000-2012 and 2012-2021 are also presented in Figure 8c,d, respectively.The elevation changes from 2000 to 2012 indicated that the whole region of the Puruogangri Ice Field was in an almost stable state, except for the thinning of several termini (e.g., those of the G94 and G98 glaciers).On the other hand, glacier thickening occurred at the terminus of the G05 Glacier.The spatial distribution of the elevation changes from 2012 to 2021 was similar to that of the elevation changes from 2000 to 2012, but the thinning magnitude was larger (Figure 9).Although thinning also occurred in the upper parts of the G13 and G98 glaciers, substantial thicknesses could be noted at their termini, which indicated that a large volume of glacial material had been transported from the upper parts to the termini of both glaciers.The mass balance of the Puruogangri Ice Field during 2000-2012 and 2012-2021 was −0.04 ± 0.01 m w.e. a −1 and −0.17 ± 0.01 m w.e. a −1 , respectively, which indicated that the glaciers were in a stable state during 2000-2012 while they were ablated significantly over the last 10 years.
In order to determine the possible stable periods and accelerated ablation periods for the Puruogangri Ice Field over the past 10 years, we further measured the thickness changes during 2012-2015 (Figure 8e), 2015-2017 (Figure 8f) and 2017-2021 (Figure 8g).Only the G94 glacier changed significantly during 2012-2015, while the rest of the ice field remained in a nearly stable state in both the accumulation and ablation regions.Similarly, this glacier was also in a stable state during 2015-2017.On the contrary, an extensive range of spatial variations were found during 2017-2021, including significant changes at the termini of several large glaciers (e.g., G94, G13 and G98).The transport of material along the G98 glacier resulted in significant thinning at the upper end and significant thickening at the lower end.Figure 9 exhibits the glacier mass balance values of −0.12 ± 0.02 m w.e. a −1 , −0.03 ± 0.01 m w.e. a −1 and −0.40 ± 0.03 m w.e. a −1 for 2012-2015, 2015-2017 and 2017-2021, respectively, which demonstrate that the significant glacier changes after 2012 were mainly concentrated in the period of 2017-2021.
In this study, we analyzed the area changes of six selected glaciers that had distinctive variation characteristics (as discussed in Section 4.1) and we calculated the corresponding results for the mass balance of these glaciers during different periods, as shown in Figure 10.The mass balance of the selected glaciers showed three types of characteristics: a continuous loss of mass (i.e., the G94, G37 and G36 glaciers, as shown in Figure 10a-c); a loss of mass and then accumulation (i.e., the G05 and G13 glaciers, as shown in Figure 10d, e); and a loss of mass and then near-stable state (i.e., the G98 glacier, as shown in Figure 10f).As examples of the first type of characteristics, the G94, G37 and G36 glaciers were ablated during the periods of both 1975-2000 and 2000-2021, for which the mass balance values were comparable to that of the G94 glacier over the same two periods, while the loss rate decreased over the two most recent decades for both the G37 and G36 glaciers, whereas the G36 glacier was in an almost stable state during the last two decades.Compared to 2000-2012, the mass balance in both the G94 and G37 glaciers decreased during 2012-2021.As examples of the second type of characteristics, the mass balance of the G05 and G13 glaciers changed from mass loss   In this study, we analyzed the area changes of six selected glaciers that had distinctive variation characteristics (as discussed in Section 4.1) and we calculated the corresponding results for the mass balance of these glaciers during different periods, as shown in Figure 10.The mass balance of the selected glaciers showed three types of characteristics: a continuous loss of mass (i.e., the G94, G37 and G36 glaciers, as shown in Figure 10ac); a loss of mass and then accumulation (i.e., the G05 and G13 glaciers, as shown in Figure 10d, e); and a loss of mass and then near-stable state (i.e., the G98 glacier, as shown in Figure 10f).As examples of the first type of characteristics, the G94, G37 and G36 glaciers were ablated during the periods of both 1975-2000 and 2000-2021, for which the mass balance values were comparable to that of the G94 glacier over the same two periods, while the loss rate decreased over the two most recent decades for both the G37 and G36 glaciers, whereas the G36 glacier was in an almost stable state during the last two decades.

Assessment of DEM Accuracy
In this study, we used ICESat-1 altimetry data to evaluate the accuracy of the DEMs that we used, based on the differences between the ICESat-1 footprints and the corresponding values of the DEMs that were extracted using bilinear interpolation.Figure 11 shows the distribution histograms of the elevation differences between the ICESat-1 data and the values of the KH-9, SRTM C-band (SRTM-C), Copernicus and ZY-3 DEMs.As can be seen from the figure, the mean elevation differences were 1.03 m, 3.24 m, −0.08 m, −2.35

Assessment of DEM Accuracy
In this study, we used ICESat-1 altimetry data to evaluate the accuracy of the DEMs that we used, based on the differences between the ICESat-1 footprints and the corresponding values of the DEMs that were extracted using bilinear interpolation.Figure 11 shows the distribution histograms of the elevation differences between the ICESat-1 data and the values of the KH-9, SRTM C-band (SRTM-C), Copernicus and ZY-3 DEMs.As can be seen from the figure, the mean elevation differences were 1.03 m, 3.24 m, −0.08 m, −2.35 m, 17.53 m and 2.78 m for the KH-9 DEM, SRTM DEM, Copernicus DEM and ZY-3 DEMs, respectively.In particular, a large average bias was found in the ZY-3 DEM that was acquired in 2017, which was mainly caused by a planar trend and was corrected in the subsequent data processing.The results showed that the accuracy of the KH-9 DEM was lower than those of the other five DEMs, as indicated by its larger standard deviation (STD).The STD of the KH-9 DEM was 8.62 m, whereas those of the SRTM-C, Copernicus and ZY-3 DEMs from 2015 were all less than 5 m.Although the accuracy of the KH-9 DEM was slightly coarse, it has been demonstrated that such an accuracy would still be suitable for quantifying glacier elevation changes over a decade timescale [26,58].

Comparison to Previous Studies
Table 3 summarizes the existing mass balance measurements for the Puruogangri Ice Field, from which it can be seen that studies that were based on RS data focused on the period of 2000-2016.Our study extended this period to 1975-2021 to help to monitor more

Comparison to Previous Studies
Table 3 summarizes the existing mass balance measurements for the Puruogangri Ice Field, from which it can be seen that studies that were based on RS data focused on the period of 2000-2016.Our study extended this period to 1975-2021 to help to monitor more recent changes in the glaciers in the Puruogangri Ice Field.For the period of 1975-2000, our mass balance result of −0.23 ± 0.02 m w.e. a −1 , which was based on satellite data, was basically consistent with that of Lei et al. [59] (approximately −0.20 m w.e. a −1 ), who mainly used topographic maps and SRTM DEMs.During the period of 2000-2012, our mass balance result of −0.04 ± 0.01 m w.e. a −1 was comparable to those of Liu et al. [30] (from −0.035 to +0.019 m w.e. a −1 , with an average of −0.019 m w.e. a −1 ), Huintjes et al. [28] (slight mass loss of −0.04 m w.e. a −1 ) and Neckel et al. [29] (−0.04 m w.e. a −1 ).Our result was slightly different from the latter two studies, which could be related to the effects of the differences between the data that were utilized and the methods that were adopted.Liu et al. [60] reported an annual surface thinning rate of −0.31 ± 0.03 m w.e. a −1 during 2012-2016 and Liu et al. [21] derived annual glacier mass balance values of 0.44 ± 0.10, −0.13 ± 0.03, −0.34 ± 0.06 and −0.52 ± 0.10 m w.e. a −1 for 2012-2012, 2012-2013, 2013-2014 and 2014-2016, respectively, while our result for 2012-2015 was −0.09 ± 0.02 m w.e. a −1 .In summary, our results were generally consistent with those of previous studies and were reliable.However, we found that a surging phenomenon occurred at the G98 glacier in 2020, as discussed in Section 4.2, which also proved that our study on recent periods was necessary.Table 3.A comparison of our study to previous studies on the glacier mass balance of the Puruogangri Ice Field.

Potential Impacts of Climate Change on Glacier Changes in the Puruogangri Ice Field
Based on the temperature and precipitation data that were recorded by the three nearest meteorological stations (Bange (31 • 13 48 N, 90 • 0 36 E), Tuotuohe (34 • 7 48 N, 92 • 15 36 E) and Anduo (32 • 12 36 N, 91 • 3 36 E) at the altitudes of 4700 m, 4533 m and 4800 m, respectively), we analyzed climate change in this region (Figure 12).From 1975 to 1999, the average annual precipitation levels that were observed at the three stations were 321.45

Conclusions
In this study, a KH-9 image and Landsat TM/ETM+/OLI images were used to estimate changes in the glacier areas of the Puruogangri Ice Field and KH-9, SRTM-C/-X, Copernicus and ZY-3 DEMs were employed to monitor the changes in glacier surface elevation and mass balance using the geodetic method.The results showed that the total area of glaciers in the Puruogangri Ice Field reduced from 427.44 ± 12.The trends of precipitation slightly increased, with annual cumulative precipitation increasing at a rate of 9.1 mm/10 a, 17.8 mm/10 a and 8.0 mm/10 a, respectively, which also showed large interannual fluctuations.We found an increase in temperature from 1975 to 2018, with the mean annual air temperature at the Bange, Tuotuohe and Anduo stations increasing at a rate of 0.51 • C/10 a, 0.55 • C/10 a and 0.36 • C/10 a, respectively.Our results for the glacier mass balance during the two periods were −0.23 ± 0.02 and −0.25 ± 0.02 m w.e. a −1 and the monitored glacier mass loss increased with the temperature.It could be inferred that the accelerated glacier mass loss in the Puruogangri Ice Field was more sensitive to rising temperatures and that our findings were consistent with those of Zhang et al. [25] for the interior of the TP.

Conclusions
In this study, a KH-9 image and Landsat TM/ETM+/OLI images were used to estimate changes in the glacier areas of the Puruogangri Ice Field and KH-9, SRTM-C/-X, Copernicus and ZY-3 DEMs were employed to monitor the changes in glacier surface elevation and mass balance using the geodetic method.The results showed that the total area of glaciers in the Puruogangri Ice Field reduced from 427.44 ± 12.43 km 2 in 1975 to 387.87 ± 11.02 km 2 in 2021, which corresponded to a rate of area loss of −0.86 km 2 a −1 .Simultaneously, a significant glacier surface thinning was observed from 1975 to 2021 and the mass loss accelerated over the past nearly five decades.Specifically, at an multidecade timescale, the glacier mass balance values were −0.23 ± 0.02 m w.e. a −1 and −0.29 ± 0.02 m w.e. a −1 for 1975-2000 and 2000-2021, respectively, which indicated a sustained and stable rate of mass loss.Nevertheless, an interesting phenomenon occurred be-tween 2000 and 2012 in that the Puruogangri Ice Field glaciers remained in an approximately steady state, with a mass balance of −0.04 ± 0.01 m w.e. a −1 .In contrast, the glaciers, overall, began to shift to a negative mass balance after 2012 and the mass balance was −0.17 ± 0.01 m w.e. a −1 during 2012-2021.Furthermore, under more intensive data sampling (i.e., 2-to 4-year timescales), we found that the region-wide glacier mass balance values were −0.12 ± 0.02 m w.e. a −1 , −0.03 ± 0.01 m w.e. a −1 and −0.40 ± 0.03 m w.e. a −1 , for 2012-2015, 2015-2017 and 2017-2021, respectively.This demonstrated that the shift in glacier state in the Puruogangri Ice Field began in 2017.It is worth pointing out that a new surging event possibly occurred in October 2020.Finally, our analysis of field-based meteorological observation data indicated that the glacier mass loss and area loss in the Puruogangri Ice Field region could be related to temperature increase.Our study provides an overview of the overall changes in the Puruogangri Ice Field and could help to reveal the responses of glaciers to long-term climate change.

Figure 1 .
Figure 1.An overview of the Puruogangri Ice Field in the central part of the TP.The background is the SRTM DEM.

Figure 1 .
Figure 1.An overview of the Puruogangri Ice Field in the central part of the TP.The background is the SRTM DEM.

24 Figure 2 .
Figure 2. The workflow of our glacier area and mass balance measurements.

Figure 2 .
Figure 2. The workflow of our glacier area and mass balance measurements.3.1.2.ZY-3 DEMs ZY-3 stereo images that were collected in 2015, 2017 and 2021 were employed to derive DEMs using the OrthoEngine module in PCI Geomatica software from NAD, FWD and BWD images.The procedures that were used to generate DEMs using the OrthoEngine software were as follows.We created a new project, input the images, collected the GCPs/TPs, computed the model, created the epipolar images and then extracted the DEMs automatically.Following the method that was used by Kropáček et al. [49], King et al. [50] and Robson et al. [51], we manually collected approximately 40 GCPs from Google Earth images and automatically collected 128 TPs for each of the three ZY-3 images.The root mean square error (RMSE) values of the GCPs and TPs were both constrained to less than 3 m.

24 Figure 3 .
Figure 3.The distribution of radar penetration depths according to altitude.The penetration depths were measured by differencing the SRTM-X and SRTM-C DEMs.3.4.Estimation of Changes in Glacier Elevation and Mass Balance 3.4.1.Estimation of Glacier Elevation Changes

Figure 3 .
Figure 3.The distribution of radar penetration depths according to altitude.The penetration depths were measured by differencing the SRTM-X and SRTM-C DEMs.

Figure 4 .
Figure 4.The distribution of normalized elevation differences in terrain aspect over stable reg (a) before and (b) after DEM co-registration.

Figure 4 .
Figure 4.The distribution of normalized elevation differences in terrain aspect over stable regions (a) before and (b) after DEM co-registration.

Figure 5 .
Figure 5.The elevation differences and the corresponding statistical histograms over stable regions (a) before and (b) after corrections for the period of 1975-2000.

Figure 5 .
Figure 5.The elevation differences and the corresponding statistical histograms over stable regions (a) before and (b) after corrections for the period of 1975-2000.

Figure 6 .
Figure 6.(a) The glacier boundaries in 1975, 2000, 2012, 2015, 2017 and 2021 (the background is the Landsat 8 image from 2021); (b) the rates of area change in the Puruogangri Ice Field over the different periods.

Figure 6 .
Figure 6.(a) The glacier boundaries in 1975, 2000, 2012, 2015, 2017 and 2021 (the background is the Landsat 8 image from 2021); (b) the rates of area change in the Puruogangri Ice Field over the different periods.

Figure 9 .
Figure 9.A comparison of the mass balance of the Puruogangri Ice Field during different periods.
Compared to 2000-2012, the mass balance in both the G94 and G37 glaciers decreased during 2012-2021.As examples of the second type of characteristics, the mass balance of the G05 and G13 glaciers changed from mass loss over 1975-2000 (−0.04 ± 0.01 m w.e. a −1 and −0.04 ± 0.01 m w.e. a −1 , respectively) to mass accumulation over 2000-2021 (0.06 ± 0.02 m w.e. a −1 and 0.03 ± 0.01 m w.e. a −1 , respectively).Compared to 2000-2012, the accumulation rate showed decreasing and increasing trends over 2012-2021 for the G05 and G13 glaciers, respectively.As an example of the third type of characteristics, the mass balance of the G98 glacier during 1975-2000 and 2000-2021 was −0.25 ± 0.02 m w.e. a −1 and −0.06 ± 0.02 m w.e. a −1 , respectively, which implied that a considerable reduction in mass loss occurred after 2000 compared to the prior period.

Figure 9 . 24 Figure 10 .
Figure 9.A comparison of the mass balance of the Puruogangri Ice Field during different periods.Remote Sens. 2022, 14, x FOR PEER REVIEW 18 of 24

Figure 12 .
Figure 12.The annual precipitation levels and average temperatures that were recorded by the meteorological stations at Bange, Tuotuohe and Anduo from 1975 to 2018.The black dashed line indicates the linear fitting change in precipitation and the black solid line indicates the linear fitting change in temperature.
43 km 2 in 1975 to 387.87 ± 11.02 km 2 in 2021, which corresponded to a rate of area loss of −0.86 km 2 a −1 .Simultaneously, a significant glacier surface thinning was observed from 1975 to 2021 and the mass loss accelerated over the past nearly five decades.Specifically, at an multidecade timescale, the glacier mass balance values were −0.23 ± 0.02 m w.e. a −1 and −0.29 ± 0.02 m w.e. a −1 for 1975-2000 and 2000-2021, respectively, which indicated a sustained and stable rate of mass loss.Nevertheless, an interesting phenomenon occurred between 2000 and 2012 in that the Puruogangri Ice Field glaciers remained in an approximately steady state, with a mass balance of −0.04 ± 0.01 m w.e. a −1 .In contrast, the glaciers, overall, began to shift to a negative mass balance after 2012 and the mass balance was −0.17 ± 0.01 m w.e. a −1 during 2012-2021.Furthermore, under more intensive data sampling (i.e., 2-to 4-year timescales), we found that the region-wide glacier mass balance values were −0.12 ± 0.02 m w.e. a −1 , −0.03 ± 0.01 m w.e. a −1 and −0.40 ± 0.03 m w.e. a −1 , for 2012-2015, 2015-2017 and 2017-2021, respectively.This demonstrated that the shift in glacier state in the Puruogangri Ice Field began in 2017.It is worth pointing out that a new surging event possibly occurred in October 2020.Finally, our analysis of field-based meteorological observation data indicated that the glacier mass loss and area loss in the Puruogangri Ice Field region could be related to temperature increase.Our study provides an overview of the overall changes in the Puruogangri Ice Field and could help to reveal the responses of glaciers to long-term cli-

Figure 12 .
Figure 12.The annual precipitation levels and average temperatures that were recorded by the meteorological stations at Bange, Tuotuohe and Anduo from 1975 to 2018.The black dashed line indicates the linear fitting change in precipitation and the black solid line indicates the linear fitting change in temperature.

Table 1 .
The data that were used in this study.

Table 1 .
The data that were used in this study.