Next Article in Journal
Using PRISMA Hyperspectral Satellite Imagery and GIS Approaches for Soil Fertility Mapping (FertiMap) in Northern Morocco
Next Article in Special Issue
Mass Balance of the Antarctic Ice Sheet in the Early 21st Century
Previous Article in Journal
Integration Data Model of the Bathymetric Monitoring System for Shallow Waterbodies Using UAV and USV Platforms
Previous Article in Special Issue
Arctic Multiyear Ice Areal Flux and Its Connection with Large-Scale Atmospheric Circulations in the Winters of 2002–2021
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
College of Earth and Environmental Sciences, Lanzhou University, Lanzhou 730000, China
2
National Tibetan Plateau Data Center (TPDC), State Key Laboratory of Tibetan Plateau Earth System, Environment and Resources (TPESER), Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100101, China
3
Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(16), 4078; https://doi.org/10.3390/rs14164078
Submission received: 18 July 2022 / Revised: 15 August 2022 / Accepted: 18 August 2022 / Published: 20 August 2022
(This article belongs to the Special Issue Remote Sensing of Ice Loss Tracking at the Poles)

Abstract

:
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 field within the Tibetan Plateau, the Puruogangri Ice Field has attracted a lot of attention from the scientific community. However, relevant studies that are based on satellite data have mainly focused on a few periods from 2000–2016. Long-term and multiperiod observations remain to be conducted. To this end, we estimated the changes in the glacier area and mass balance of the Puruogangri Ice Field over five subperiods between 1975 and 2021, based on multisource remote sensing data. Specifically, we employed KH-9 and Landsat images to estimate the area change from 1975 to 2021 using the band ratio method. Subsequently, based on KH-9 DEM, SRTM DEM, Copernicus DEM and ZY-3 DEM data, we evaluated the glacier elevation changes and mass balance over five subperiods during 1975–2021. The results showed that the total glacier area decreased from 427.44 ± 12.43 km2 to 387.87 ± 11.02 km2 from 1975 to 2021, with a decrease rate of 0.86 km2 a−1. The rate of area change at a decade timescale was −0.74 km2 a−1 (2000–2012) and −1.00 km2 a−1 (2012–2021). Furthermore, the rates at a multiyear timescale were −1.23 km2 a−1, −1.83 km2 a−1 and −0.42 km2 a−1 for 2012–2015, 2015–2017 and 2017–2021, respectively. In terms of the glacier mass balance, the region-wide results at a two-decade timescale were −0.23 ± 0.02 m w.e. a−1 for 1975–2000 and −0.29 ± 0.02 m w.e. a−1 for 2000–2021, indicating a sustained and relatively stable mass loss over the past nearly five decades. After 2000, the loss rate at a decade timescale was −0.04 ± 0.01 m w.e. a−1 for 2000–2012 and −0.17 ± 0.01 m w.e. a−1 for 2012–2021, indicating an increasing loss rate over recent decades. It was further found that the mass loss rate was −0.12 ± 0.02 m w.e. a−1 for 2012–2015, −0.03 ± 0.01 m w.e. a−1 for 2015–2017 and −0.40 ± 0.03 m w.e. a−1 for 2017–2021. These results indicated that a significant portion of the glacier mass loss mainly occurred after 2017. According to our analysis of the meteorological measurements in nearby regions, the trends of precipitation and the average annual air temperature both increased. Combining these findings with the results of the glacier changes implied that the glacier changes seemed to be more sensitive to temperature increase in this region. Overall, our results improved our understanding of the status of glacier changes and their reaction to climate change in the central Tibetan Plateau.

1. 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 × 106 km2 and the volume of total ice reserves is about 27.6 × 106 km3, 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 km2 (7.37%) from 1992 to 2009 and shrank by 19.90 km2 (2.84 km2 a−1), 39.95 km2 (5.71 km2 a−1) and 6.83 km2 (2.28 km2 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 km2 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.

2. Study Area and Data

2.1. 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 km2, 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 northwest–southeast 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].

2.2. Data

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

2.2.1. Optical Satellite Images

(1)
Landsat
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 long-term 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 km2 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).

2.2.2. 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).

2.2.3. 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).

2.2.4. Meteorological Data

Since there is no weather station on the Puruogangri Ice Field, the temperature and precipitation records were collected from the three adjacent meteorological stations in 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) (Figure 1 and Table 1). The daily recorded temperature and precipitation data for the period of 1975–2018 were collected in this study, which are available at https://data.cma.cn/ (last access: 11 March 2021).

3. 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.

3.1. DEM Extraction

3.1.1. 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.

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.

3.2. 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.

3.3. 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.

3.4. Estimation of Changes in Glacier Elevation and Mass Balance

3.4.1. 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:
h 1 = a · cos ( b β ) · tan α + d h ¯  
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 ( d h ¯ ) 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.
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)):
h 2 = a 1 x 2 + a 2 y 2 + a 3 x y + a 4 x + a 5 y + a 6
h 3 = b 1 z 2 + b 2 z + b 3
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.

3.4.2. 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:
B = i = 1 n ( h i ¯ · A i ) · ρ g A · ρ w
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).

3.5. Uncertainty Assessment

The uncertainty of the glacier boundary delineation was estimated using the length of the glacier boundaries and the image resolution [10]:
δ = 1 2 · l · r
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:
δ h = δ E N e f f
N e f f = N t · R 2 D
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 ( δ m 0 ) 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):
δ m 0 = ( S · ρ g ρ w · S t · δ h ) 2 + ( h · ρ g ρ w · S t · δ s ) 2 + ( S · h ρ w · S t · δ ρ g ) 2
δ m = δ m 0 2 + δ C X 2 + δ X 2
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.

4. Results

4.1. 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 km2 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 km2 a−1, with rates of −0.83 km2 a−1 and −0.90 km2 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 km2 a−1 for 2000–2012 and −1.31 km2 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 km2 a−1, −2.74 km2 a−1 and 0.40 km2 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 Figure 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 km2 (−0.05 km2 a−1), −3.83 km2 (−0.15 km2 a−1) and −0.73 km2 (−0.03 km2 a−1) between 1975 and 2000, respectively, and −2.20 km2 (−0.10 km2 a−1), −3.04 km2 (−0.14 km2 a−1) and −1.07 km2 (−0.05 km2 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 km2 during 1975–2000, expanded 0.53 km2 during 2000–2015 and then started to shrink again, decreasing from 14.62 ± 0.31 km2 (2015) to 14.40 ± 0.31 km2 (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 km2 a−1 and −0.14 km2 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 km2 (+0.34 km2 a−1). The area of the G13 Glacier also slightly increased by +0.32 km2 (+0.08 km2 a−1) during this period.

4.2. 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 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.

5. Discussion

5.1. 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].

5.2. 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.

5.3. 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 mm, 267.40 mm and 436.37 mm, respectively, while from 2000 to 2018, they were 348.05 mm, 329.56 mm and 462.97 mm, respectively. The average temperature was −0.75, −4.15 and −2.50 °C from 1975 to 1999 and 0.43, −2.67 and −1.56 °C from 2000 to 2018, respectively.
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.

6. 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 km2 in 1975 to 387.87 ± 11.02 km2 in 2021, which corresponded to a rate of area loss of −0.86 km2 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 climate change.

Author Contributions

Conceptualization, S.R., X.L. and D.Z.; methodology, S.R., Y.Z. and Y.W.; software, S.R., Y.Z. and Y.W.; validation, S.R., Y.Z. and Y.W.; formal analysis, S.R., Y.Z. and D.Z.; investigation, X.L., D.Z. and D.J.; resources, X.L., D.Z. and D.J.; data curation, S.R.; writing—original draft preparation, S.R.; writing—review and editing, D.Z., Y.Z. and Y.W.; visualization, S.R., Y.Z. and Y.W.; supervision, X.L., D.Z. and Y.N.; project administration, X.L., D.Z. and Y.N.; funding acquisition, Y.Z., X.L. and D.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the State Key Research and Development Program of China (grant number 2021YFB3900105) and the National Natural Science Foundation of China (grant numbers 41988101 and 42001381).

Data Availability Statement

Not applicable.

Acknowledgments

We would like to thank the USGS (https://earthexplorer.usgs.gov/, accessed on 26 October 2021, 19 August 2022 and 2 October 2021) for providing the Landsat images and the SRTM-C DEM and KH-9 stereo images. The ZY-3 stereo images were from the Natural Resources Satellite Remote Sensing Cloud Service Platform (http://sasclouds.com/chinese/normal/, accessed on 23 January 2022) and the SRTM-X DEM was from DLR (https://eoweb.dlr.de/egp/, accessed on 19 August 2022). The ICESat-1 altimetry data were from the US National Snow and Ice Data Center (https://nsidc.org/data/GLAH14/versions/34, accessed on 23 October 2021). The daily climate data were from the China Meteorological Data Service Center (https://data.cma.cn/, last access on 11 March 2021).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bolch, T.; Kulkarni, A.; Kääb, A.; Huggel, C.; Paul, F.; Cogley, J.G.; Frey, H.; Kargel, J.S.; Fujita, K.; Scheel, M.; et al. The state and fate of himalayan glaciers. Science 2012, 336, 310–314. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Shean, D.E.; Bhushan, S.; Montesano, P.; Rounce, D.R.; Arendt, A.; Osmanoglu, B. A systematic, regional assessment of high mountain asia glacier mass balance. Front. Earth Sci. 2020, 7, 363. [Google Scholar] [CrossRef] [Green Version]
  3. Cao, B.; Guan, W.; Li, K.; Wen, Z.; Han, H.; Pan, B. Area and mass changes of glaciers in the west kunlun mountains based on the analysis of multi-temporal remote sensing images and dems from 1970 to 2018. Remote Sens. 2020, 12, 2632. [Google Scholar] [CrossRef]
  4. Yang, J.; Li, M.; Yang, S.; Tan, C. Vulnerability of the glaciers to climate change in china:Current situation and evaluation. J. Glaciol. Geocryol. 2013, 35, 1077–1087. [Google Scholar]
  5. Yao, T.; Wu, G.; Xu, B.; Wang, W.; Gao, J.; An, B. Asian water tower change and its impacts. Bull. Chin. Acad. Sci. 2019, 34, 1203–1209. [Google Scholar]
  6. Gao, J.; Yao, T.; Masson-Delmotte, V.; Steen-Larsen, H.C.; Wang, W. Collapsing Glaciers Threaten Asia’s Water Supplies; Nature Publishing Group: Beijing, China, 2019. [Google Scholar]
  7. Jacob, T.; Wahr, J.; Pfeffer, W.T.; Swenson, S. Recent contributions of glaciers and ice caps to sea level rise. Nature 2012, 482, 514–518. [Google Scholar] [CrossRef]
  8. Zhou, Y.; Li, X.; Zheng, D.; Li, Z.; An, B.; Wang, Y.; Jiang, D.; Su, J.; Cao, B. The joint driving effects of climate and weather changes caused the chamoli glacier-rock avalanche in the high altitudes of the india himalaya. Sci. China Earth Sci. 2021, 64, 1909–1921. [Google Scholar] [CrossRef]
  9. Kääb, A.; Leinss, S.; Gilbert, A.; Bühler, Y.; Gascoin, S.; Evans, S.G.; Bartelt, P.; Berthier, E.; Brun, F.; Chao, W.-A. Massive collapse of two glaciers in western tibet in 2016 after surge-like instability. Nat. Geosci. 2018, 11, 114–120. [Google Scholar] [CrossRef] [Green Version]
  10. Wang, Y.; Li, J.; Wu, L.; Guo, L.; Hu, J.; Zhang, X. Estimating the changes in glaciers and glacial lakes in the xixabangma massif, central himalayas, between 1974 and 2018 from multisource remote sensing data. Remote Sens. 2021, 13, 3903. [Google Scholar] [CrossRef]
  11. Veh, G.; Korup, O.; von Specht, S.; Roessner, S.; Walz, A. Unchanged frequency of moraine-dammed glacial lake outburst floods in the himalaya. Nat. Clim. Change 2019, 9, 379–383. [Google Scholar] [CrossRef]
  12. Zhou, Y.; Hu, J.; Li, Z.; Li, J.; Zhao, R.; Ding, X. Quantifying glacier mass change and its contribution to lake growths in central kunlun during 2000–2015 from multi-source remote sensing data. J. Hydrol. 2019, 570, 38–50. [Google Scholar] [CrossRef]
  13. Zhou, Y.; Li, Z.; Li, J.I.A. Slight glacier mass loss in the karakoram region during the 1970s to 2000 revealed by kh-9 images and srtm dem. J. Glaciol. 2017, 63, 331–342. [Google Scholar] [CrossRef] [Green Version]
  14. Zhou, Y.; Li, Z.; Li, J.I.A.; Zhao, R.; Ding, X. Geodetic glacier mass balance (1975–1999) in the central pamir using the srtm dem and kh-9 imagery. J. Glaciol. 2019, 65, 309–320. [Google Scholar] [CrossRef] [Green Version]
  15. Wang, N.; Yao, T.; Xu, B.; Chen, A.A.; Wang, W. Spatiotemporal pattern, trend, and influence of glacier change in tibetan plateau and surroundings under global warming. Bull. Chin. Acad. Sci. 2019, 34, 1220–1232. [Google Scholar]
  16. Zhang, J.; He, X.; Ye, B.; Wu, J. Recent variation of mass balance of the xiao dongkemadi glacier in the tanggula range and its influencing factors. J. Glaciol. Geocryol. 2013, 35, 263–271. [Google Scholar]
  17. The surface flow features of the puruogangri ice field. J. Glaciol. Geocryol. 2003, 25, 288–290.
  18. Kumar, A.; Negi, H.S.; Kumar, K.; Shekhar, C.; Kanda, N. Quantifying mass balance of east-karakoram glaciers using geodetic technique. Polar Sci. 2019, 19, 24–39. [Google Scholar] [CrossRef]
  19. Gardelle, J.; Berthier, E.; Arnaud, Y.; Kääb, A. Region-wide glacier mass balances over the pamir-karakoram-himalaya during 1999–2011. Cryosphere 2013, 7, 1263–1286. [Google Scholar] [CrossRef] [Green Version]
  20. Hugonnet, R.; McNabb, R.; Berthier, E.; Menounos, B.; Nuth, C.; Girod, L.; Farinotti, D.; Huss, M.; Dussaillant, I.; Brun, F.; et al. Accelerated global glacier mass loss in the early twenty-first century. Nature 2021, 592, 726–731. [Google Scholar] [CrossRef]
  21. Liu, L.; Jiang, L.; Jiang, H.; Wang, H.; Ma, N.; Xu, H. Accelerated glacier mass loss (2011–2016) over the puruogangri ice field in the inner tibetan plateau revealed by bistatic insar measurements. Remote Sens. Environ. 2019, 231, 111241. [Google Scholar] [CrossRef]
  22. Zhou, W.; Li, Z.; Li, J.; Cui, Z.; Wang, C. Variations of glaciers and glacial lake in geladandong mountain range in 1992–2009 with remote-sensing technology. J. Cent. South Univ. 2014, 2014, 45. [Google Scholar]
  23. Ye, Q.; Kang, S.; Chen, F.; Wang, J. Monitoring glacier variations on geladandong mountain, central tibetan plateau, from 1969 to 2002 using remote-sensing and gis technologies. J. Glaciol. 2006, 52, 537–545. [Google Scholar] [CrossRef] [Green Version]
  24. Berthier, E.; Brun, F. Karakoram geodetic glacier mass balances between 2008 and 2016: Persistence of the anomaly and influence of a large rock avalanche on siachen glacier. J. Glaciol. 2019, 65, 494–507. [Google Scholar] [CrossRef] [Green Version]
  25. Zhang, Z.; Jiang, L.; Liu, L.; Sun, Y.; Wang, H. Annual glacier-wide mass balance (2000–2016) of the interior tibetan plateau reconstructed from modis albedo products. Remote Sens. 2018, 10, 1031. [Google Scholar] [CrossRef] [Green Version]
  26. Zhou, Y.; Li, Z.; Li, J.; Zhao, R.; Ding, X. Glacier mass balance in the qinghai–tibet plateau and its surroundings from the mid-1970s to 2000 based on hexagon kh-9 and srtm dems. Remote Sens. Environ. 2018, 210, 96–112. [Google Scholar] [CrossRef]
  27. Liu, G.; Fan, J.-h.; Zhao, F.; Mao, K.-b. Monitoring elevation change of glaciers on geladandong mountain using tandem-x sar interferometry. J. Mt. Sci. 2017, 14, 859–869. [Google Scholar] [CrossRef]
  28. Huintjes, E.; Neckel, N.; Hochschild, V.; Schneider, C. Surface energy and mass balance at purogangri ice cap, central tibetan plateau, 2001–2011. J. Glaciol. 2015, 61, 1048–1060. [Google Scholar] [CrossRef] [Green Version]
  29. Neckel, N.; Braun, A.; Kropáček, J.; Hochschild, V. Recent mass balance of the purogangri ice cap, central tibetan plateau, by means of differential x-band sar interferometry. Cryosphere 2013, 7, 1623–1633. [Google Scholar] [CrossRef] [Green Version]
  30. Liu, L.; Jiang, L.; Sun, Y.; Wang, H.; Yi, C.; Hsu, H. Morphometric controls on glacier mass balance of the puruogangri ice field, central tibetan plateau. Water 2016, 8, 496. [Google Scholar] [CrossRef] [Green Version]
  31. Thompson, L.G.; Tandong, Y.; Davis, M.E.; Mosley-Thompson, E.; Mashiotta, T.A.; Lin, P.-N.; Mikhalenko, V.N.; Zagorodnov, V.S. Holocene climate variability archived in the puruogangri ice cap on the central tibetan plateau. Ann. Glaciol. 2006, 43, 61–69. [Google Scholar] [CrossRef] [Green Version]
  32. Maurer, J.; Rupper, S. Tapping into the hexagon spy imagery database: A new automated pipeline for geomorphic change detection. ISPRS J. Photogramm. Remote Sens. 2015, 108, 113–127. [Google Scholar] [CrossRef]
  33. Surazakov, A.; Aizen, V. Positional accuracy evaluation of declassified hexagon kh-9 mapping camera imagery. Photogramm. Eng. Remote Sens. 2010, 76, 603–608. [Google Scholar] [CrossRef] [Green Version]
  34. Dong, Y.; Chen, W.; Chang, H.; Zhang, Y.; Feng, R.; Meng, L. Assessment of orthoimage and dem derived from zy-3 stereo image in northeastern china. Surv. Rev. 2016, 48, 247–257. [Google Scholar] [CrossRef]
  35. Hasegawa, H.; Matsuo, K.; Koarai, M.; Watanabe, N.; Fukushima, Y. Dem accuracy and the base to height (b/h) ratio of stereo images. Int. Arch. Photogramm. Remote Sens. 2000, 33, 356–359. [Google Scholar]
  36. Farr, T.G.; Rosen, P.A.; Caro, E.; Crippen, R.; Duren, R.; Hensley, S.; Kobrick, M.; Paller, M.; Rodriguez, E.; Roth, L.; et al. The shuttle radar topography mission. Rev. Geophys. 2007, 45. [Google Scholar] [CrossRef] [Green Version]
  37. Rabus, B.; Eineder, M.; Roth, A.; Bamler, R. The shuttle radar topography mission—A new class of digital elevation models acquired by spaceborne radar. ISPRS J. Photogramm. Remote Sens. 2003, 57, 241–262. [Google Scholar] [CrossRef]
  38. AIRBUS. Copernicus Dem: Copernicus Digital Elevation Model Product Handbook; AIRBUS: Leiden, The Netherlands, 2020. [Google Scholar]
  39. Karlson, M.; Bastviken, D.; Reese, H. Error characteristics of pan-arctic digital elevation models and elevation derivatives in northern sweden. Remote Sens. 2021, 13, 4653. [Google Scholar] [CrossRef]
  40. Purinton, B.; Bookhagen, B. Beyond vertical point accuracy: Assessing inter-pixel consistency in 30 m global dems for the arid central andes. Front. Earth Sci. 2021, 9, 758606. [Google Scholar] [CrossRef]
  41. Brosens, L.; Campforts, B.; Govers, G.; Aldana-Jague, E.; Razanamahandry, V.F.; Razafimbelo, T.; Rafolisy, T.; Jacobs, L. Comparative analysis of the copernicus, tandem-x, and uav-sfm digital elevation models to estimate lavaka (gully) volumes and mobilization rates in the lake alaotra region (madagascar). Earth Surf. Dynam. 2022, 10, 209–227. [Google Scholar] [CrossRef]
  42. Round, V.; Leinss, S.; Huss, M.; Haemmig, C.; Hajnsek, I. Surge dynamics and lake outbursts of kyagar glacier, karakoram. Cryosphere 2017, 11, 723–739. [Google Scholar] [CrossRef] [Green Version]
  43. Lambrecht, A.; Mayer, C.; Wendt, A.; Floricioiu, D.; VÖLksen, C. Elevation change of fedchenko glacier, pamir mountains, from gnss field measurements and tandem-x elevation models, with a focus on the upper glacier. J. Glaciol. 2018, 64, 637–648. [Google Scholar] [CrossRef] [Green Version]
  44. Marešová, J.; Gdulová, K.; Pracná, P.; Moravec, D.; Gábor, L.; Prošek, J.; Barták, V.; Moudrý, V. Applicability of data acquisition characteristics to the identification of local artefacts in global digital elevation models: Comparison of the copernicus and tandem-x dems. Remote Sens. 2021, 13, 3931. [Google Scholar] [CrossRef]
  45. Yan, L.; Wang, J.; Shao, D. Glacier mass balance in the manas river using ascending and descending pass of sentinel 1a/1b data and srtm dem. Remote Sens. 2022, 14, 1506. [Google Scholar] [CrossRef]
  46. Zhao, F.; Long, D.; Li, X.; Huang, Q.; Han, P. Rapid glacier mass loss in the southeastern tibetan plateau since the year 2000 from satellite observations. Remote Sens. Environ. 2022, 270, 112853. [Google Scholar] [CrossRef]
  47. Kropáček, J.; Neckel, N.; Bauder, A. Estimation of mass balance of the grosser aletschgletscher, swiss alps, from icesat laser altimetry data and digital elevation models. Remote Sens. 2014, 6, 5614–5632. [Google Scholar] [CrossRef] [Green Version]
  48. Sochor, L.; Seehaus, T.; Braun, M.H. Increased ice thinning over svalbard measured by icesat/icesat-2 laser altimetry. Remote Sens. 2021, 13, 2089. [Google Scholar] [CrossRef]
  49. Kropáček, J.; Vařilová, Z.; Baroň, I.; Bhattacharya, A.; Eberle, J.; Hochschild, V. Remote sensing for characterisation and kinematic analysis of large slope failures: Debre sina landslide, main ethiopian rift escarpment. Remote Sens. 2015, 7, 16183–16203. [Google Scholar] [CrossRef] [Green Version]
  50. King, O.; Quincey, D.J.; Carrivick, J.L.; Rowan, A.V. Spatial variability in mass loss of glaciers in the everest region, central himalayas, between 2000 and 2015. Cryosphere 2017, 11, 407–426. [Google Scholar] [CrossRef] [Green Version]
  51. Robson, B.A.; Nuth, C.; Nielsen, P.R.; Girod, L.; Hendrickx, M.; Dahl, S.O. Spatial variability in patterns of glacier change across the manaslu range, central himalaya. Front. Earth Sci. 2018, 6, 12. [Google Scholar] [CrossRef] [Green Version]
  52. Pieczonka, T.; Bolch, T.; Junfeng, W.; Shiyin, L. Heterogeneous mass loss of glaciers in the aksu-tarim catchment (central tien shan) revealed by 1976 kh-9 hexagon and 2009 spot-5 stereo imagery. Remote Sens. Environ. 2013, 130, 233–244. [Google Scholar] [CrossRef] [Green Version]
  53. Zhou, Y.; Li, X.; Zheng, D.; Zhang, X.; Wang, Y.; Ren, S.; Guo, Y. Decadal changes in glacier area, surface elevation and mass balance for 2000–2020 in the eastern tanggula mountains using optical images and tandem-x radar data. Remote Sens. 2022, 14, 506. [Google Scholar] [CrossRef]
  54. Li, C.; Jiang, L.; Liu, L.; Wang, H. Regional and altitude-dependent estimate of the srtm c/x-band radar penetration difference on high mountain asia glaciers. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2021, 14, 4244–4253. [Google Scholar] [CrossRef]
  55. Nuth, C.; Kääb, A. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. Cryosphere 2011, 5, 271–290. [Google Scholar] [CrossRef] [Green Version]
  56. Bolch, T.; Pieczonka, T.; Benn, D.I. Multi-decadal mass loss of glaciers in the everest area (nepal himalaya) derived from stereo imagery. Cryosphere 2011, 5, 349–358. [Google Scholar] [CrossRef] [Green Version]
  57. Huss, M. Density assumptions for converting geodetic glacier volume change to mass change. Cryosphere 2013, 7, 877–887. [Google Scholar] [CrossRef] [Green Version]
  58. Zhang, X.; Zhou, J.; Liu, Z. Dem extraction and precision evaluation of mountain glaciers in the qianghai-tibet plateau based on kh-9 data: Take the purog kangri glacier and the jiong glacier as example. J. Glaciol. Geocryol. 2019, 41, 27–35. [Google Scholar]
  59. Lei, Y.; Yao, T.; Yi, C.; Wang, W.; Sheng, Y.; Li, J.; Joswiak, D. Glacier mass loss induced the rapid growth of linggo co on the central tibetan plateau. J. Glaciol. 2012, 58, 177–184. [Google Scholar] [CrossRef] [Green Version]
  60. Liu, L.; Jiang, L.; Sun, Y.; Yi, C.; Wang, H.; Hsu, H. Glacier elevation changes (2012–2016) of the puruogangri ice field on the tibetan plateau derived from bi-temporal tandem-x insar data. Int. J. Remote Sens. 2016, 37, 5687–5707. [Google Scholar] [CrossRef]
Figure 1. An overview of the Puruogangri Ice Field in the central part of the TP. The background is the SRTM DEM.
Figure 1. An overview of the Puruogangri Ice Field in the central part of the TP. The background is the SRTM DEM.
Remotesensing 14 04078 g001
Figure 2. The workflow of our glacier area and mass balance measurements.
Figure 2. The workflow of our glacier area and mass balance measurements.
Remotesensing 14 04078 g002
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 3. The distribution of radar penetration depths according to altitude. The penetration depths were measured by differencing the SRTM-X and SRTM-C DEMs.
Remotesensing 14 04078 g003
Figure 4. The distribution of normalized elevation differences in terrain aspect over stable regions (a) before and (b) after DEM co-registration.
Figure 4. The distribution of normalized elevation differences in terrain aspect over stable regions (a) before and (b) after DEM co-registration.
Remotesensing 14 04078 g004
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. The elevation differences and the corresponding statistical histograms over stable regions (a) before and (b) after corrections for the period of 1975–2000.
Remotesensing 14 04078 g005
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. (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.
Remotesensing 14 04078 g006
Figure 7. A regional display of the area changes in six typical glaciers: (a) the G94 glacier; (b) the G37 glacier; (c) the G36 glacier; (d) the G05 glacier; (e) the G13 glacier; (f) the G98 glacier.
Figure 7. A regional display of the area changes in six typical glaciers: (a) the G94 glacier; (b) the G37 glacier; (c) the G36 glacier; (d) the G05 glacier; (e) the G13 glacier; (f) the G98 glacier.
Remotesensing 14 04078 g007
Figure 8. Maps of the glacier surface elevation changes: (a) 1975–2000; (b) 2000–2021; (c) 2000–2012; (d) 2012–2021; (e) 2012–2015; (f) 2015–2017; (g) 2017–2021; the mass balance values are presented in the last panel (bottom right).
Figure 8. Maps of the glacier surface elevation changes: (a) 1975–2000; (b) 2000–2021; (c) 2000–2012; (d) 2012–2021; (e) 2012–2015; (f) 2015–2017; (g) 2017–2021; the mass balance values are presented in the last panel (bottom right).
Remotesensing 14 04078 g008
Figure 9. A comparison of the mass balance of the Puruogangri Ice Field during different periods.
Figure 9. A comparison of the mass balance of the Puruogangri Ice Field during different periods.
Remotesensing 14 04078 g009
Figure 10. Comparisons of the mass balance of six selected glaciers over different periods: (a) the G94 glacier; (b) the G36 glacier; (c) the G37 glacier; (d) the G05 glacier; (e) the G13 glacier; (f) the G98 glacier.
Figure 10. Comparisons of the mass balance of six selected glaciers over different periods: (a) the G94 glacier; (b) the G36 glacier; (c) the G37 glacier; (d) the G05 glacier; (e) the G13 glacier; (f) the G98 glacier.
Remotesensing 14 04078 g010
Figure 11. The distribution histograms of the elevation differences between the ICESat-1 data and the (a) KH-9, (b) SRTM-C, (c) Copernicus and (d) ZY-3 DEMs from 2015, (e) the ZY-3 DEM from 2017 and (f) the ZY-3 DEM from 2021.
Figure 11. The distribution histograms of the elevation differences between the ICESat-1 data and the (a) KH-9, (b) SRTM-C, (c) Copernicus and (d) ZY-3 DEMs from 2015, (e) the ZY-3 DEM from 2017 and (f) the ZY-3 DEM from 2021.
Remotesensing 14 04078 g011
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.
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.
Remotesensing 14 04078 g012
Table 1. The data that were used in this study.
Table 1. The data that were used in this study.
DataAcquisition DateResolutionPurpose
Landsat 5/TM7 November 200030 mGlacier boundary delineation for 2000
Landsat 7/ETM+29 September 2012Pan: 15 mGlacier boundary delineation for 2012
Landsat 8/OLI30 September 2015
21 October 2017
17 November 2021
Pan: 15 mGlacier boundary delineation for 2015, 2017 and 2021
KH-9 Image20 December 19756–9 mGlacier boundary delineation and DEM extraction for 1975
ZY-3 Stereo Images30 December 2015
15 November 2017
24 November 2021
NAD: 2.1 m; DLC: 3.5 mDEM extraction for 2015, 2017 and 2021
SRTM DEMs
(C- and X-bands)
~ February
2000
30 mDEM acquisition for 2000 and estimation of penetration depth for C-band radar
Copernicus DSM~July 201230 mDEM acquisition for 2012
ICESat-1/GLAH142003–2009~70 mDEM accuracy assessment
Meteorological
Data
1975–2018Station Climate change analysis
Table 2. The glacier areas in 1975, 2000, 2012, 2015, 2017 and 2021.
Table 2. The glacier areas in 1975, 2000, 2012, 2015, 2017 and 2021.
Glacier
(GLIMS_ID)
Glacier Area (km2)
197520002012201520172021
Puruogangri Ice Field427.44 ± 12.43406.78 ± 11.65399.70 ± 11.38392.56 ± 11.11387.08 ± 11.20387.87 ± 11.02
G9450.13 ± 0.6548.87 ± 0.6247.65 ± 0.6247.34 ± 0.6246.86 ± 0.6146.67 ± 0.61
G3743.86 ± 0.8140.03 ± 0.7537.62 ± 0.7337.42 ± 0.7737.07 ± 0.7636.99 ± 0.75
G3613.11 ± 0.3012.38 ± 0.2911.54 ± 0.2711.49 ± 0.2611.35 ± 0.2611.31 ± 0.25
G0514.28 ± 0.2914.09 ± 0.2914.50 ± 0.3114.62 ± 0.3114.42 ± 0.3114.40 ± 0.31
G139.06 ± 0.248.57 ± 0.238.42 ± 0.238.40 ± 0.198.27 ± 0.238.59 ± 0.24
G9830.43 ± 0.6729.93 ± 0.6125.90 ± 0.6025.25 ± 0.5424.57 ± 0.6425.91 ± 0.63
Table 3. A comparison of our study to previous studies on the glacier mass balance of the Puruogangri Ice Field.
Table 3. A comparison of our study to previous studies on the glacier mass balance of the Puruogangri Ice Field.
ResearchDataPeriodResult
Lei et al. [59]Topographic maps and SRTM DEMs1974–2000−0.20 m w.e. a−1 a
Huintjes et al. [28]High Asia Refined analysis (HAR) dataset2000–2011−0.038 m w.e. a−1 b
Neckel et al. [29]SRTM-X and TanDEM-X DEMs2000–2012−0. 044 ± 0.002 m w.e. a−1
Liu et al. [30]SRTM-X and TanDEM-X DEMs2000–2012−0.035 to +0.019 m w.e. a−1
Liu et al. [60]TanDEM-X DEMs2012–2016−0.317 ± 0.027 m w.e. a−1
Liu et al. [21]TanDEM-X DEMs2011–20160.44 ± 0.10 m w.e. a−1 over 2011–2012
−0.13 ± 0.03 m w.e. a−1 over 2012–2013
−0.34 ± 0.06 m w.e. a−1 over 2013–2014
−0.52 ± 0.10 m w.e. a−1 over 2014–2016
a, The elevation change result that was presented by the authors was −6.19 m, which was converted to −0.20 m w.e. a−1 using the mass balance calculation formula; b, the annual mass balance result that was presented by the authors was −44 kg m−2 a−1, which was converted to −0.044 m w.e. a−1 by dividing by the water density (1000 kg m−3).
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Ren, S.; Li, X.; Wang, Y.; Zheng, D.; Jiang, D.; Nian, Y.; Zhou, Y. Multitemporal Glacier Mass Balance and Area Changes in the Puruogangri Ice Field during 1975–2021 Based on Multisource Satellite Observations. Remote Sens. 2022, 14, 4078. https://doi.org/10.3390/rs14164078

AMA Style

Ren S, Li X, Wang Y, Zheng D, Jiang D, Nian Y, Zhou Y. Multitemporal Glacier Mass Balance and Area Changes in the Puruogangri Ice Field during 1975–2021 Based on Multisource Satellite Observations. Remote Sensing. 2022; 14(16):4078. https://doi.org/10.3390/rs14164078

Chicago/Turabian Style

Ren, Shanshan, Xin Li, Yingzheng Wang, Donghai Zheng, Decai Jiang, Yanyun Nian, and Yushan Zhou. 2022. "Multitemporal Glacier Mass Balance and Area Changes in the Puruogangri Ice Field during 1975–2021 Based on Multisource Satellite Observations" Remote Sensing 14, no. 16: 4078. https://doi.org/10.3390/rs14164078

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop