Estimating the Changes in Glaciers and Glacial Lakes in the Xixabangma Massif, Central Himalayas, between 1974 and 2018 from Multisource Remote Sensing Data

: The continuous melting of valley glaciers can impact the water levels of glacial lakes and create glacial lake outburst ﬂoods (GLOFs). The Xixabangma massif is one of the most populated areas in the Himalayas and has suffered from multiple GLOFs. To estimate the glacier melting rate in the past four decades and analyze the outburst risk of glacial lakes in the Xixabangma massif, we determined changes in glacier mass balance, glacier area and glacial lake area based on KH-9 images, TanDEM-X images, Landsat images, SRTM DEM and ICESat-2 elevations. Our results show that, from 1974 to 2018, the total glacier area shrank from 954.01 km 2 to 752.46 km 2 , whereas the total glacial lake area grew from 20.90 km 2 to 38.71 km 2 . From 1974 to 2000, 2000 to 2013 and 2013 to 2018, the region-wide glacier mass balance values were − 0.16 m w.e./a, − 0.31 m w.e./a and − 0.29 m w.e./a, respectively. Three glacial lakes, named Gangxico, Galongco and Jialongco, respectively, expanded by 127.14%, 373.45% and 436.36% from 1974 to 2018, and the mass loss rates of their parent glaciers from 2000 to 2013 increased by 81.72%, 122.22% and 160.00% relative to those during 1974 to 2000. The dams of these three lakes are unstable, and their drainage valleys directly connect to a major town and its infrastructure. Due to current high-water levels, possible external events such as ice collapse, landslide, heavy rainfall and earthquakes can easily trigger GLOFs. Hence, we deemed that the Gangxico, Galongco and Jialongco glacial lakes are dangerous and require special attention.


Introduction
The Himalayas host more than 12,000 glaciers that cover an area of approximately 33,000 km 2 and store approximately 12,000 km 3 of fresh water [1,2]. Abundant glacial meltwater enables the Himalayas to form large rivers, such as the Brahmaputra, the Ganges and the Indus [3]. Since the 1950s, climate warming has aggravated glacier ablation in the Himalayas [4], and many proglacial and supraglacial lakes, which are located in the upstream parts of valleys, have emerged after glaciers have experienced constant backwasting and retreating. As temperature maintains its high level, glacier mass loss will continue and more meltwater will gather in glacial lakes. Worse still, some glaciers connected to proglacial lakes will experience an accelerated retreat, leaving more space for glacial lake expansion. Consequently, water pressure over lake dams will increase, leading to the potential hazard of dam bursts. Hence, glacier mass loss can directly affect the stability of a glacial lake, so the quantitative estimation of glacier mass balance can help to assess the hazard of glacial lake outbursts [5]. At present, glacier mass balance can be measured in several ways, including in situ measurements (e.g., stake records, differential global In order to reveal the evolution of glaciers and glacial lakes in the Xixabangma massif during the past 40 years, and to analyze the risk of glacial lake outbursts based on multiple factors including glacier mass balance, we first delineated four terms related to glacier and Remote Sens. 2021, 13, 3903 3 of 21 glacial lake boundaries and derived geodetic glacier mass balance in three periods during 1974-2018 from multi-source remote sensing data. Then, we analyzed glaciers' responses to climate change between 1967 and 2018. Finally, the risk of glacial lake outbursts was discussed based on the expanding rates and sizes of lakes, the mass balance of parent glaciers and other current factors. plentiful precipitation, and the glaciers developed there are primarily maritime in nature. Since the transport of water vapor is blocked by high mountains, the northern slope receives much less precipitation, and the glaciers developed there belong primarily to the continental type. The annual average temperature is approximately 2.1 °C, and the annual accumulated precipitation is approximately 580 mm [27]. Precipitation mainly occurs in summer. Many large glaciers have their tongues covered by debris, and rivers are fed by glacier meltwater flow in the region's winding and deep-cut valleys. Proglacial moraine and outwash are abundant, and GLOFs can easily create debris flow. The resident population of the Xixabangma massif is approximately 25,000. Many infrastructures, including the China-Nepal Road, have been built in this region. This region is close to the Himalayan seismic belt, which suffers from frequent earthquakes; an earthquake up to a magnitude of 7.8 (M 7.8) occurred on 25 April 2015, and another one up to M 7.3 occurred on 12 May 2015. Both triggered considerable debris flow in China and Nepal [28].

Datasets and Methods
We first delineated glacier outlines and extracted DEMs from remote sensing images (Table 2) and then derived changes in glacier thickness by differentiating the DEMs of different dates and removing systematic bias in the elevation difference maps. Afterward, we computed changes in glacier mass based on glacier area, elevation change and mass

Study Area
The Xixabangma massif (85 • 15 -86 • 05 E, 28 • 05-28 • 45 N) straddles the China-Nepal border ( Figure 1). Its altitude ranges from approximately 1000 m a.s.l. to 8000 m a.s.l. Influenced by the Indian monsoon and westerly circulation, the southern slope receives plentiful precipitation, and the glaciers developed there are primarily maritime in nature. Since the transport of water vapor is blocked by high mountains, the northern slope receives much less precipitation, and the glaciers developed there belong primarily to the continental type. The annual average temperature is approximately 2.1 • C, and the annual accumulated precipitation is approximately 580 mm [27]. Precipitation mainly occurs in summer. Many large glaciers have their tongues covered by debris, and rivers are fed by glacier meltwater flow in the region's winding and deep-cut valleys. Proglacial moraine and outwash are abundant, and GLOFs can easily create debris flow. The resident population of the Xixabangma massif is approximately 25,000. Many infrastructures, including the China-Nepal Road, have been built in this region. This region is close to the Himalayan seismic belt, which suffers from frequent earthquakes; an earthquake up to a magnitude of 7.8 (M 7.8) occurred on 25 April 2015, and another one up to M 7.3 occurred on 12 May 2015. Both triggered considerable debris flow in China and Nepal [28].

Datasets and Methods
We first delineated glacier outlines and extracted DEMs from remote sensing images (Table 2) and then derived changes in glacier thickness by differentiating the DEMs of different dates and removing systematic bias in the elevation difference maps. Afterward, we computed changes in glacier mass based on glacier area, elevation change and mass density. The workflow of glacier mass balance measurement is illustrated in Figure 2. Detailed processing is introduced in the next sections. density. The workflow of glacier mass balance measurement is illustrated in Figure 2. Detailed processing is introduced in the next sections.

Generation of KH-9 DEM
Reconnaissance optical satellite Keyhole (KH-9) was on active service from 1971 to 1986. The archive images were declassified in 2002. The KH-9 images distributed by the United States Geological Survey were scanning images of the film negative taken by a frame-mapping camera. The metadata of KH-9 imaging are unknown. The quality of released KH-9 images is affected by pixel distortion and scratches due to their long storage time. Nevertheless, given their high spatial resolution (6-9 m), wide swath (120 × 250 km),

Generation of KH-9 DEM
Reconnaissance optical satellite Keyhole (KH-9) was on active service from 1971 to 1986. The archive images were declassified in 2002. The KH-9 images distributed by the United States Geological Survey were scanning images of the film negative taken by a framemapping camera. The metadata of KH-9 imaging are unknown. The quality of released KH-9 images is affected by pixel distortion and scratches due to their long storage time. Nevertheless, given their high spatial resolution (6-9 m), wide swath (120 × 250 km), high forward overlapping rate (70%) and the advantage of early acquisition time, these KH-9 images are an important data source for land-cover change detection [29]. Each KH-9 image consists of two sub-images, and each sub-image has a reseau cross alignment of 23 × 24. In this study, one KH-9 stereo image of 7 microns (3600 dpi) acquired on 23 November 1974 was used to generate a DEM (Table 2). To extract a DEM from a KH-9 stereo image, three steps, including pre-processing, internal orientation and external orientation, need to be performed. Image pre-processing consists of image resampling, reseau cross identification, distortion correction, sub-image mosaic and contrast enhancement [30]. Internal orientation was accomplished on the basis of the geometric model of the frame-mapping camera and parameter settings. External orientation was accomplished on the basis of external control points and tie points. The external control points were manually selected at the places of ground features, such as river confluences and giant rocks, by combining Landsat-7/ETM+ L1T images and SRTM-C DEM (in ice-free regions) [31]. The tie points were automatically generated. Finally, a KH-9 DEM with a resolution of 30 m was generated.

Generation of TanDEM-X DEM
Given their high spatial resolution, wide swath and all-weather working capability, active synthetic aperture radar (SAR) images are widely used for ground deformation monitoring and topography mapping. Shuttle Radar Topography Mission (SRTM), a famous SAR interferometry (InSAR) mission that works in a single-pass dual-antenna mode, obtained DEMs covering 80% of the global land surface within 11 days in February 2000. The German Aerospace Center launched the TerraSAR-X and TanDEM-X satellites in 2007 and 2010, respectively. These two satellites form a constellation. Under the bistatic working mode, one of the two satellites transmits signals, and two receive an echo simultaneously; therefore, master and slave images for interferometry can be obtained at the same time. The zero-temporal baseline ensures the temporal correlation of the interferometric phase. In this study, three bistatic strip map TanDEM-X image pairs acquired in 2013 and 2018 were used to generate DEMs ( Table 2).
The raw interferometric phase (ϕ int ) comprises the topographic phase (ϕ topo ), the flat earth phase (ϕ f lat ), the displacement phase (ϕ dis ), the atmospheric phase (ϕ atm ) and the noise phase (ϕ noise ) Equation (1). Given that master and slave images are acquired simultaneously, the displacement and atmospheric phase can be ignored. Most of the noise phase can be removed by multilooking and filtering, and the flat earth phase can be simulated on the basis of imaging parameters. After the removal of the noise phase and the flat earth phase, the topographic phase (ϕ topo ) can be estimated.
Given the short wavelength of the X-band, the fringes of TanDEM-X interferograms are dense in the mountainous region, thus making phase unwrapping difficult. To address this problem, a topographic phase was simulated from SRTM-C DEM and subtracted from the observed TanDEM-X topographic phase. The differential phase with much sparser fringes was easily unwrapped. The unwrapped phase was converted into elevation difference through the following formula: where λ , R, θ 0 and B ⊥ are the radar wavelength, slant range, nominal incident angle and perpendicular baseline, respectively. A new DEM was obtained by adding elevation difference to the SRTM-C DEM. After backward geocoding, we finally derived the TanDEM-X DEMs with a spatial resolution of 30 m in 2013 and 2018.

Delineation of Glacier and Glacial Lake Outlines
Glacier outlines are a prerequisite for calculating glacier area. Moreover, when evaluating DEM accuracy and correcting systematic bias in elevation difference maps, we need to mask the glacierized area with glacier outlines. We collected Landsat L1T images with little cloud and seasonal snow cover and conducted radiation calibration and atmospheric correction on them [32]. The initial glacier outline (Figure 3a) was extracted from the corrected Landsat images through the band ratio method (Red/SWIR) [33]. Water bodies and fresh snowpacks are generally misclassified as glaciers, and glaciers in the shadow cannot be discerned. We manually edited the glacier outlines by combining multiple images acquired in the same seasons. Regarding debris cover in the glacier ablation zone, we further manually revised the glacier boundary based on Landsat panchromatic-multispectral fused images, Google Earth high-resolution images, DEM and elevation difference maps [23].
where , , and are the radar wavelength, slant range, nominal incident angle and perpendicular baseline, respectively. A new DEM was obtained by adding elevation difference to the SRTM-C DEM. After backward geocoding, we finally derived the Tan-DEM-X DEMs with a spatial resolution of 30 m in 2013 and 2018.

Delineation of Glacier and Glacial Lake Outlines
Glacier outlines are a prerequisite for calculating glacier area. Moreover, when evaluating DEM accuracy and correcting systematic bias in elevation difference maps, we need to mask the glacierized area with glacier outlines. We collected Landsat L1T images with little cloud and seasonal snow cover and conducted radiation calibration and atmospheric correction on them [32]. The initial glacier outline (Figure 3a) was extracted from the corrected Landsat images through the band ratio method (Red/SWIR) [33]. Water bodies and fresh snowpacks are generally misclassified as glaciers, and glaciers in the shadow cannot be discerned. We manually edited the glacier outlines by combining multiple images acquired in the same seasons. Regarding debris cover in the glacier ablation zone, we further manually revised the glacier boundary based on Landsat panchromatic-multispectral fused images, Google Earth high-resolution images, DEM and elevation difference maps [23].
The orthorectified KH-9 images were used to extract glacier outlines in 1974. We first delineated the glacier outlines from KH-9 images through visual interpretation and then manually revised them on the basis of the elevation difference map (SRTM-C DEM -KH-9 DEM).
The boundaries of glacial lakes were first extracted through a normalized difference water index method and then manually edited through visual interpretation.  The orthorectified KH-9 images were used to extract glacier outlines in 1974. We first delineated the glacier outlines from KH-9 images through visual interpretation and then manually revised them on the basis of the elevation difference map (SRTM-C DEM-KH-9 DEM).
The boundaries of glacial lakes were first extracted through a normalized difference water index method and then manually edited through visual interpretation.

Evaluation of DEM Accuracy
We used ICESat-2/ATLAS ATL06 elevation products as a reference to evaluate the accuracy of DEMs generated in this study. The accuracy of the ICESat-1/GLAS elevation product in flat regions reaches 0.14 m [34], and the officially announced accuracy of the ICESat-2/ATLAS elevation product is higher [35]. The footprints of ICESat-2 data adopted Remote Sens. 2021, 13, 3903 7 of 21 in this study are shown in Figure 4f. Regarding the negative effect of steep slopes, we selected the ICESat-2 points in stable regions with a slope below 30 degrees. The DEM values at the footprints of ICESat-2 data were extracted via bilinear interpolation. Given that outliers may exist in ICESat-2 elevation products and DEMs, we set a threshold of ±100 m for elevation difference. As shown in Figure 4, the accuracy of KH-9 DEM was lower than that of the SRTM-C DEM and TanDEM-X DEM. The root mean square error of KH-9 is ±21.53 m, whereas that for TanDEM-X and SRTM-C DEMs ranges from only ±3 m to ±7 m. However, glacier mass balance changed greatly during the last decades of the 20th century, and the differencing of KH-9 DEM (1974) with SRTM-C (2000) or TanDEM-X (2018) DEMs can still provide valuable information regarding changes in glacier thickness.

Evaluation of DEM Accuracy
We used ICESat-2/ATLAS ATL06 elevation products as a reference to evaluate the accuracy of DEMs generated in this study. The accuracy of the ICESat-1/GLAS elevation product in flat regions reaches 0.14 m [34], and the officially announced accuracy of the ICESat-2/ATLAS elevation product is higher [35]. The footprints of ICESat-2 data adopted in this study are shown in Figure 4f. Regarding the negative effect of steep slopes, we selected the ICESat-2 points in stable regions with a slope below 30 degrees. The DEM values at the footprints of ICESat-2 data were extracted via bilinear interpolation. Given that outliers may exist in ICESat-2 elevation products and DEMs, we set a threshold of ±100 m for elevation difference. As shown in Figure 4, the accuracy of KH-9 DEM was lower than that of the SRTM-C DEM and TanDEM-X DEM. The root mean square error of KH-9 is ±21.53 m, whereas that for TanDEM-X and SRTM-C DEMs ranges from only ±3 m to ±7 m. However, glacier mass balance changed greatly during the last decades of the 20th century, and the differencing of KH-9 DEM (1974) with SRTM-C (2000) or TanDEM-X (2018) DEMs can still provide valuable information regarding changes in glacier thickness.

DEM Co-Registration and Systematic Error Correction
Changes in glacier thickness were estimated by differencing DEMs at different times. Due to the different acquisition times and imaging approaches, and the different methods of DEM generation, the processed KH-9, SRTM-C and TanDEM-X DEMs have different geolocation accuracies and may cover different areas, even though they have the same geographical coordinates. Any two DEMs must be co-registered before differencing them because a tiny shift in the DEMs can cause considerable errors in elevation differences in mountainous regions. We used the analytical co-registration method proposed by Nuth and Kääb in 2011 [36]. This method is operated on the basis of terrain slope, aspect and raw elevation difference and has been proven to be robust and efficient in matching DEMs over glacierized regions. This method decomposes the DEM shift vector into vertical and horizontal components. Vertical shift is assumed to be systematic. A functional relationship is built amongst horizontal shift, terrain slope, aspect and raw elevation difference Equation (3). The shift parameters were solved on the basis of the least square principle, and the DEMs were resampled in accordance with the shift parameters. For additional details about this method, readers can refer to [36]. The distribution of normalized

DEM Co-Registration and Systematic Error Correction
Changes in glacier thickness were estimated by differencing DEMs at different times. Due to the different acquisition times and imaging approaches, and the different methods of DEM generation, the processed KH-9, SRTM-C and TanDEM-X DEMs have different geolocation accuracies and may cover different areas, even though they have the same geographical coordinates. Any two DEMs must be co-registered before differencing them because a tiny shift in the DEMs can cause considerable errors in elevation differences in mountainous regions. We used the analytical co-registration method proposed by Nuth and Kääb in 2011 [36]. This method is operated on the basis of terrain slope, aspect and raw elevation difference and has been proven to be robust and efficient in matching DEMs over glacierized regions. This method decomposes the DEM shift vector into vertical and horizontal components. Vertical shift is assumed to be systematic. A functional relationship is built amongst horizontal shift, terrain slope, aspect and raw elevation difference Equation (3). The shift parameters were solved on the basis of the least square principle, and the DEMs were resampled in accordance with the shift parameters. For additional details about this method, readers can refer to [36]. The distribution of normalized elevation difference (TanDEM-X DEM in 2013-RTM-C DEM in 2000) against terrain aspect before DEM co-registration is illustrated in Figure 5a. where dh is the elevation difference in stable region; a is the plane shift vector; b is the direction of plane shift vector; ϕ is the terrain aspect; α is the terrain slope; Z is the vertical bias.

Correction of C-band and X-band Penetration Effects
The penetration of radar signals over glacier surfaces must be considered when estimating changes in glacier thickness through differencing InSAR and optical DEMs or differencing InSAR DEMs of distinct microwave bands. In theory, the absolute X-band penetration depth can be estimated by differencing the TanDEM-X elevations with the ICTSat-2/ATLAS ATL06 elevations because the penetration depth of laser signals over glacier surfaces is negligible. However, ICEsat-2 data were not available until October 2018, whereas our newest TanDEM-X data were acquired in February 2018. We assumed the elevation change between 24 February 2018 and 24 February 2019 was equal to the average annual change between 24 February 2013 and 24 February 2018; we then extrapolated the TanDEM-X glacier surface elevation on 24 February 2019. The ICESat-2/ATLAS ATL06 elevations on 27 December 2018 and 9 February 2019 were then differenced with the extrapolated TanDEM-X elevation to estimate the winter X-band penetration depth over glacier surfaces. The distinction between C-band and X-band penetration depths over glacier surfaces was estimated by differencing the simultaneously obtained C-band Elevation difference maps usually contain systematic biases that can considerably affect the result of glacier thickness change measurements [36,37]. We also found systematic biases in the elevation difference maps that were related to planimetric position, altitude and terrain maximum curvature. In particular, planimetric position-related bias may be caused by baseline residue (for InSAR DEM generation) and spatial interferometric decorrelation (for InSAR DEM generation) [36]; altitude related bias may be caused by the non-uniform distribution of external control points in altitudes (for optical DEM generation) [38] and spatial interferometric decorrelation (for InSAR DEM generation); and terrain maximum curvature bias may be caused by a resolution discrepancy between the source data of two DEMs [37,39]. Assuming the systematic biases of elevation differences over glaciers and surrounding ice-free regions are common, we fitted these systematic biases as universal trends from elevation differences in ice-free regions and removed them from the whole elevation difference maps [36]. Systematic bias related to planimetric position can be corrected by quadratic or cubic surface fitting Equation (4), and systematic biases related to altitude and terrain curvature can be corrected by quadratic or cubic polynomials Equation (5). The order n depends on the distribution of bias. The magnitude of biases related to planimetric position, altitude and curvature (TanDEM-X DEM in 2013-SRTM-C DEM in 2000) is illustrated in Figure 5b, 5c and 5d, respectively. dh = a 0 + a 1 x + a 2 y + a 3 x 2 + a 4 xy + a 5 y 2 (4) dh = a o + a 1 p + a 2 p 2 + . . . + a n p n (5) where dh is the elevation difference observation; a i is the model coefficient; x and y are the row and column number of the pixel; p is the terrain parameters (altitude or terrain curvature); n is the order of polynomial. Elevation differences over stable regions are supposed to be zero; therefore, a standard deviation (STD) of observations over stable regions can describe the quality of elevation difference maps. However, STD is sensitive to outliers. We used normalized median absolute deviation (NMAD) instead, as it is less sensitive to outliers and can be considered a robust estimator of STD [31,40]. The NMAD is computed by Equation (6). Table 3 shows the NMAD and mean of elevation difference observations in stable regions before and after DEM co-registration and systematic bias corrections. The smaller the NMAD, the higher the accuracy of elevation difference observations.
where ∆h is the elevation difference observation and n is the number of elevation difference observations.

Correction of C-Band and X-Band Penetration Effects
The penetration of radar signals over glacier surfaces must be considered when estimating changes in glacier thickness through differencing InSAR and optical DEMs or differencing InSAR DEMs of distinct microwave bands. In theory, the absolute X-band penetration depth can be estimated by differencing the TanDEM-X elevations with the ICTSat-2/ATLAS ATL06 elevations because the penetration depth of laser signals over glacier surfaces is negligible. However, ICEsat-2 data were not available until October 2018, whereas our newest TanDEM-X data were acquired in February 2018. We assumed the elevation change between 24 February 2018 and 24 February 2019 was equal to the average annual change between 24 February 2013 and 24 February 2018; we then extrapolated the TanDEM-X glacier surface elevation on 24 February 2019. The ICESat-2/ATLAS ATL06 elevations on 27 December 2018 and 9 February 2019 were then differenced with the extrapolated TanDEM-X elevation to estimate the winter X-band penetration depth over glacier surfaces. The distinction between C-band and X-band penetration depths over glacier surfaces was estimated by differencing the simultaneously obtained C-band and X-band SRTM DEM. The absolute winter C-band penetration depth was then obtained by adding the above two results together. The SRTM-X DEM does not cover the glaciers in the Xixabangma massif. Regarding the similarity of climate condition, we selected the Lugula region (red area in the inset panel in Figure 1), which is adjacent to the Xixabangma massif, to estimate the difference of C-band and X-band penetration depths. Finally, we fitted the penetration depth as a function of altitude and slope and corrected the penetration effects in elevation difference maps pixel by pixel.

Computation of Glacier Mass Balance
The change in glacier volume was computed through the hypsography method, and a mass density of 850 ± 60 kg/m 3 was used to convert the glacier volume change into glacier mass balance [41]. The computation formula of glacier mass balance is where B, A, ρ glacier and ρ water are glacier mass balance, glacier area, glacier mass density and water density (1000 kg/m 3 ), respectively; N is the number of altitude intervals; A i and ∆h i are the glacier area and average glacier thickness change within each altitude interval (100 m), respectively.

Uncertainty Analysis
Assuming the glacier/lake boundary passed through the center of each pixel, we computed the uncertainty of glacier/lake (δ area ) based on the image spatial resolution and the length of the glacier/lake boundary [42]. The computation formula is where l is the length of the glacier/lake boundary and r is the image resolution. The uncertainty of the average glacier elevation (δ average ) was calculated on the basis of the uncertainty of the elevation difference over stable regions and the spatial autocorrection of the elevation difference [7,43]. The computation formula is where N MAD is calculated by Formula (6), A is the observed glacier area and d is the spatial autocorrection distance of elevation difference in the stable area. In this study, we assumed that the spatial correlation was 30 pixels, namely 900 m. For glacier mass density, we adopted ±60 kg/m 3 as the uncertainty. The uncertainty of X-band penetration depth was ±0.48 m, and the uncertainty of the difference of X-band and C-band penetration depths was ±0.17 m. The uncertainty of glacier mass balance was determined through the basic error propagation law. The computation formula is where δ mass , δ area , δ average and δ ρ glacier are the uncertainties of glacier mass balance, glacier area, glacier elevation difference and glacier mass density, respectively; ∆h is the glacier elevation difference; A is the glacier area; S total is the total glacier area. Table 4 shows the areas of glaciers and glacial lakes in 1974, 2000, 2012 and 2018. From 1974 to 2018, the glacier area shrank by 201.55 km 2 (−0.48%/a), whereas the glacial lake area expanded by 17.81 km 2 (1.94%/a). In three consecutive periods, the glacier shrinking rate during the period from 2012 to 2018 was the highest (5.69 km 2 /a), and the glacial lake expanding rate from 2000 to 2012 was the highest (0.53 km 2 /a). The discrepancy in these changes may be caused by two factors. First, many glacial lakes have reached their maximum extent and are bounded by an abrupt steeping of topography [17]. Second, even though the glacier shrinking rate during the period from 2012 to 2018 was the highest, the glacier mass loss rate during the period from 2000 to 2012 was the highest. We primarily studied six large glacial lakes in the Xixabangma massif. The locations of the six glacial lakes are shown in Figure 6. The glacial lakes Guoqiangco and Guoruco were mainly fed by the glaciers Guoluo and Guoluoqiang, respectively. From 1974 to 2018, the glacial lakes Guoqiangco and Guoruco expanded by 0.41 km 2 (8.32%) and 0.61 km 2 (14.66%), respectively. Meanwhile, their parent glaciers shrank by 1.49 km 2 (17.37%) and 1.51 km 2 (10.55%), respectively. Glacial lakes Gangxico, Galongco and Jialongco also have parent glaciers. However, relative to glacial lakes Guoqiangco and Guoruco, these three glacial lakes expanded much faster. From 1974 to 2018, the areas of Gangxico, Galongco and Jialongco increased by 2.53 km 2 (127.41%), 4.22 km 2 (373.45%) and 0.48 km 2 (436.36%), respectively, and the area of their parent glaciers decreased by 3.23 km 2 (−56.27%), 5.01 km 2 (−25.32%) and 1.35 km 2 (−23.85%), respectively. After the year 2000, the glacial lake Galongco expanded at a rate of 0.12 km 2 /a, and its parent glacier, Jicongpu, shrank at a rate of 0.14 km 2 /a. Both rates are considerable. Glacial lake Kungco does not have a parent glacier. From 1974 to 2018, its area decreased from 2.22 km 2 to 1.92 km 2 , and its shrinking rates after 2000 were even higher than that before 2000. Table 4. Changes in glacier and glacial lake areas in the Xixabangma massif from 1974 to 2018.

Glacier (GLIMS ID)
Glacial Lake  Figure 7a shows the measured X-band penetration depths in different altitude ranges. Over ice-free areas, the X-band penetration depths are close to zero, whereas over glacierized areas, the X-band penetration depths range from 0.15 m to 5.52 m. The mean values over ice-free and glacierized areas are 0.10 and 3.79 m, respectively. Dehecq et al. [44] measured an average X-band penetration depth of 4 m at the glacier accumulation zone of the Mont-Blanc area by differencing the Pléiades DEM and TanDEM-X DEM, whereas Lambrecht et al. [45] measured an X-band penetration depth of 3-4 m at the accumulation zone of the Fedchenko Glacier, Pamir, by differentiating between the GNSS elevation and TanDEM-X elevation.  whereas over glacierized areas, the penetration depth differences range from −0.88 m to 3.16 m. The mean values over ice-free and glacierized areas are −0.25 and 1.44 m, respectively. Gardelle et al. [22] measured a penetration depth difference of 1.4 m over the western Nepal glaciers by differencing SRTM-X and SRTM-C DEM, whereas Zhou et al. [46] measured a penetration depth difference of 1.4 m over the central Nepal glaciers through the same method. The radar penetration depth over glaciers generally in-creases with altitude because colder temperature is generally more suitable for pre-serving fresh snow.

Changes in Glacier Thickness and Mass Balance
The measured region-wide glacier mass balance in the Xixabangma massif was −0.   Figure 7b shows the measured differences between X-band and C-band penetration depths. Similarly, over ice-free areas, the penetration depth differences are close to zero, whereas over glacierized areas, the penetration depth differences range from −0.88 m to 3.16 m. The mean values over ice-free and glacierized areas are −0.25 and 1.44 m, respectively. Gardelle et al. [22] measured a penetration depth difference of 1.4 m over the western Nepal glaciers by differencing SRTM-X and SRTM-C DEM, whereas Zhou et al. [46] measured a penetration depth difference of 1.4 m over the central Nepal glaciers through the same method. The radar penetration depth over glaciers generally in-creases with altitude because colder temperature is generally more suitable for pre-serving fresh snow.

Changes in Glacier Thickness and Mass Balance
The measured region-wide glacier mass balance in the Xixabangma massif was −0.  G085812E28335N) had the highest receding rate, and its terminus was directly connected to glacial lake Galongco. The glacial lake expanded upward as the glacier terminus retreated. Thus, the lengths of thickness change profiles during three periods vary (Profile C). The cumulative thinning of glacier Jicongpu during the first two consecutive periods reached 55.69 m and 41.06 m, and the corresponding annual thinning rates along the profile were 1.08 m/a and 0.99 m/a. However, the annual thickness change rate along the profile during 2013-2018 was +0.11 m/a. As shown in Figure 6, the ablation zone of glacier Jicongpu (GLIMS ID: G085812E28335N) almost vanished before 2012; hence, the observed thickness change rate during 2013-2018 represents a mass change in its accumulation zone. Moreover, glacier Jicongpu had the highest thinning rates. Table 5 lists the average thickness change rates of seven glaciers, including the largest one in the Xixabangma massif and the parent glaciers of the six glacial lakes mentioned above. change rate along the profile during 2013-2018 was +0.11 m/a. As shown in Figure 6, the ablation zone of glacier Jicongpu (GLIMS ID: G085812E28335N) almost vanished before 2012; hence, the observed thickness change rate during 2013-2018 represents a mass change in its accumulation zone. Moreover, glacier Jicongpu had the highest thinning rates. Table 5 lists the average thickness change rates of seven glaciers, including the largest one in the Xixabangma massif and the parent glaciers of the six glacial lakes mentioned above.    −0.14 ± 0.12 −0.27 ± 0.11 --

Comparison with Previous Glacier Mass Balance Measurements
To corroborate the reliability of our glacier mass balance results, we compared them with previously published results that were derived during the same periods. Before our study, Raina et al. [47]  for glaciers in the Bhutan-China border based on KH-9 and ASTER stereo images. These glacier mass balance results were not affected by radar penetration over glacier surfaces. After removing the effect of C-band penetration, we derived a glacier mass balance of −0.16 m w.e./a from 1974 to 2000, which agreed with above three results. In addition, Pellicciotti et al. [25] derived a mass balance of −0.32 m w.e./a from 1974 to 2000 for glaciers in the Nepal Langtang based on KH-9 stereo images and SRTM-C DEM, and Wu et al. [48] derived a mass balance of −0.24 m w.e./a from 1980 to 2000 for the Kangri Karpo glacier in the Everest region based on topographic maps and SRTM-C DEM. These two studies only removed the effect of relative radar penetration depth, i.e., the difference between X-band and C-band penetration depths. If we used the same method to deal with radar penetration, then a glacier mass balance of −0.26 m w.e./a from 1974 to 2000 could be derived; this value is also close to these two results. Furthermore, Kääb et al. [49]

Relationship between Climate Change and Glacier/Glacial Lake Dynamics
Relative to glacier area and length, glacier mass balance generally responds to climate change more directly. Analyzing the precipitation and temperature data recorded by the Nyalam meteorological station (black triangle in Figure 1), we found that, from 1967 to 2018, the annual average temperature in the Xixabangma massif increased at a rate of 0.23 • C/10a, whereas the annual cumulative precipitation decreased at a rate of 11.4 mm/10a. The climate in the Xixabangma massif showed a warming and drying trend over the past half century ( Figure 10).

Outburst Risk of Major Glacial Lakes in the Xixabangma Massif
Here, we discuss the outburst risk of six glacial lakes (the specific location is shown in Figure 6) in the Xixabangma massif based on the expanding rates and sizes of lakes, the mass balance of parent glaciers, the state of dams and the distribution of downstream towns.
Glacial lake Kungco is located on the east side of the Xixabangma massif. Its area is 1.92 km 2 in 2018. During the three consecutive area observation periods (1974-2000, 2000-2012 and 2012-2018), it shrank by −1.67%, −6.41% and −6.12%, respectively. It has no parent glacier, and its shrinking may be related to a lack of precipitation. Moreover, it was formed in a nature depression and has no dam. Hence, we deemed that the outburst risk of glacial lake Kungco to be low.
Glacial lakes Guoqiangco and Guoruco were large in size in 2018 (5.34 km 2 and 4.77 km 2 ). From 1974 to 2018, they expanded slowly, and their parent glaciers lost mass moderately. In particular, during the three consecutive area observation periods, Guoqiangco expanded by 2.84%, 2.76% and 2.50%, respectively, and Guoruco expanded by 12.26%, 1.50% and 0.63 %, respectively. During the first two mass balance observation periods, the mass balance of Guoqiangco's parent glacier, Guoluo, was −0.12 m w.e./a and −0.23 m w.e./a, respectively, and the corresponding values of Guoruco' parent glacier, Guoluoqiang, were −0.04 m w.e./a and −0.10 m w.e./a, respectively (the number of observations during 2013 to 2018 was too small). Currently, the dams of these two glacial lakes are about 300 m wide and formed by solid sediments. These two glacial lakes are located on the north side of the Xixabangma massif and are connected to lake Peikuco by the Naqu river. Towns and infrastructure are absent in the downstream regions. Their impacts on the safety of the residents and social properties are small. Hence, in spite of their large size, we deemed the outburst risk of glacial lakes Guoqiangco and Guoruco to be low.
Glacial lakes Gangxico and Galongco were also large in size in 2018 (4.52 km 2 and 5.35 km 2 ). During 1974 to 2018, they expanded quickly, and their parent glaciers lost mass remarkably. In particular, during the three consecutive area observation periods, Gangxico expanded by 75.69%, 28.36% and 0.67%, respectively, and Galongco expanded by A rise in temperature intensifies glacier ablation, and a decrease in precipitation reduces mass accumulation. The annual average temperatures during the three consecutive glacier mass balance observation periods (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) and 2013-2018) were 3.57 • C, 4.27 • C and 4.22 • C, and the corresponding annual cumulative precipitation was 666.5 mm, 612.2 mm and 602.8 mm, respectively. As mentioned above, our glacier mass balance results in the three consecutive periods were −0.16 m w.e./a, −0.31 m w.e./a and −0.29 m w.e./a, respectively. In this case, the changing pattern in glacier mass balance aligned with the changing pattern in temperature. This point agrees with that in our previous study [37], i.e., glacier mass balance is more sensitive to changes in temperature than to changes in precipitation.
In this study, the area of lake Kungco that has no water supply from glaciers decreased during the observation period, which is in accordance with the change in precipitation in the Xixabangma massif. By contrast, the area of lakes Guoqiangco, Guoruco, Gangxico, Jialongco and Galongco that have direct water supply from parent glaciers increased. The continuous rise in temperature in the Xixabangma massif may accelerate glacier mass melting and led to an increase in glacial lake water area. Hence, glacier mass balance should be considered a factor in the hazard assessment of glacial lakes [5].

Outburst Risk of Major Glacial Lakes in the Xixabangma Massif
Here, we discuss the outburst risk of six glacial lakes (the specific location is shown in Figure 6) in the Xixabangma massif based on the expanding rates and sizes of lakes, the mass balance of parent glaciers, the state of dams and the distribution of downstream towns.
Glacial lake Kungco is located on the east side of the Xixabangma massif. Its area is 1.92 km 2 in 2018. During the three consecutive area observation periods (1974-2000, 2000-2012 and 2012-2018), it shrank by −1.67%, −6.41% and −6.12%, respectively. It has no parent glacier, and its shrinking may be related to a lack of precipitation. Moreover, it was formed in a nature depression and has no dam. Hence, we deemed that the outburst risk of glacial lake Kungco to be low.
Glacial lakes Guoqiangco and Guoruco were large in size in 2018 (5.34 km 2 and 4.77 km 2 ). From 1974 to 2018, they expanded slowly, and their parent glaciers lost mass moderately. In particular, during the three consecutive area observation periods, Guoqiangco expanded by 2.84%, 2.76% and 2.50%, respectively, and Guoruco expanded by 12.26%, 1.50% and 0.63 %, respectively. During the first two mass balance observation periods, the mass balance of Guoqiangco's parent glacier, Guoluo, was −0.12 m w.e./a and −0.23 m w.e./a, respectively, and the corresponding values of Guoruco' parent glacier, Guoluoqiang, were −0.04 m w.e./a and −0.10 m w.e./a, respectively (the number of observations during 2013 to 2018 was too small). Currently, the dams of these two glacial lakes are about 300 m wide and formed by solid sediments. These two glacial lakes are located on the north side of the Xixabangma massif and are connected to lake Peikuco by the Naqu river. Towns and infrastructure are absent in the downstream regions. Their impacts on the safety of the residents and social properties are small. Hence, in spite of their large size, we deemed the outburst risk of glacial lakes Guoqiangco and Guoruco to be low.
Glacial lakes Gangxico and Galongco were also large in size in 2018 (4.52 km 2 and 5.35 km 2 ). During 1974 to 2018, they expanded quickly, and their parent glaciers lost mass remarkably. In particular, during the three consecutive area observation periods, Gangxico expanded by 75.69%, 28.36% and 0.67%, respectively, and Galongco expanded by 187.07%, 55.88% and 5.98%, respectively. During the first two mass balance observation periods, the mass balance of Gangxico's parent glacier, Reqiang, was −0.93 m w.e./a and −1.69 m w.e./a., respectively (no observation during 2013 to 2018), and the corresponding values of Galongco's parent glacier, Jicongpu, were −0.36 m w.e./a and −0.80 m w.e./a, respectively. Note that the mass balance of Jicongpu during the third period seemed to be abnormal because its ablation zone almost melted away at the beginning of this period. Currently, the dams of these two glacial lakes are about 400 m and 300 m wide, respectively. Field investigations show that their dams are composed of stagnant ice and moraine [16]. A rise in temperature could cause the ice to melt inside the dam, thus weakening its stability. The downstream drainage channels of these two lakes are directly connected to Nyalam County and the China-Nepal road, one of the main traffic roads between two countries. Moreover, substantial loose debris is found in the downstream drainage channel. The GLOFs can easily create debris flow and cause extensive damage to downstream infrastructures and residential areas. Hence, we deemed that the outburst risk of glacial lakes Gangxico and Galongco to be high.
Glacial lake Jialongco was small in size in 2018 (0.59 km 2 ). Between 1974 and 2018, it expanded quickly, and its parent glacier, Lengbugang, lost mass remarkably. In particular, during the three consecutive area observation periods, its area increased by 76.58%, 190.82% and 3.51%, respectively. During the three consecutive mass balance observation periods, the mass balance values of Lengbugang were −0.10 m w.e./a, −0.26 m w.e./a and −0.28 m w.e./a, respectively. This glacial lake is located in a high place and the outside slope of its dam is steep. Currently, its dam is about 50 m wide. Water has already overflowed at the lowest part of the dam and discharge flux has scoured out a shallow channel. The drainage channel of Jialongco is connected to Nyalam County and the distance is much shorter than those found in Gangxico and Galongco. Similarly, there is also much loose debris in the downstream drainage channel, and a lake outburst flood could easily turn into debris flow. Due to strong glacier melting and heavy rainfall, this lake burst twice in 2002, causing severe losses in Nyalam County and on the China-Nepal international road. If temperatures keep rising, we expect more glacial meltwater to enter the lake, and the channel in the dam will enlarge as scouring intensifies; in turn, this may cause the dam to burst again. Hence, we deemed that the outburst risk of lake Jialongco to be very high.
In 1950s and 1960s, when the glacier melting had not been intensified yet, GLOFs also happened in the Himalayas, because external events such as ice collapse, landslide, and heavy rainfall can directly trigger the lake outburst. However, the increase in lake water volume caused by the serious glacier mass loss during the recent 20 years aggravated the situation undoubtedly. As the volume of lake water increases, the stability of the lake dam become increasingly susceptible to the external events mentioned above. Besides, as shown in Figure 6, glaciers Lengbugang, Reqiang and Jicongpu significantly retreated after 2000, especially for the latter two, whose terminus directly connect with the glacial lakes. Without buttressing from the lower reaches, the higher reaches of the glaciers will collapse more easily. Furthermore, the Xixabangma massif is located in seismic active zone. Earthquake not only can cause ice collapse and landslides, but also can damage the lake dams directly [28]. In general, glacial lakes Gangxico, Galongco and Jialongco are dangerous and needs special attention.

Limitations of This Study
The water volume of a glacial lake is determined by its area and depth. However, we have no information on lake depth. Given the weak applicability of current area-depth scaling approaches [51], lake depth should be measured for the evaluation of future glacial lake risks. Meanwhile, we analyzed the outburst risk of glacial lakes in a qualitative manner. Compared with quantitative evaluation methods, such as fuzzy consistent matrix [52], event tree [53] and additive ratio scales [54], a qualitative analysis based on multiple factors is somewhat subjective. More rigorous risk evaluation results should be obtained through the use of quantitative methods.
Besides, we have no field data other than meteorological data from one climate station. Given the topographic variation, the data collected by the climate station far away from the glaciers may not accurately characterize the local temperature and precipitation. The glacier debris cover and debris flux also play important roles on regulating glacier ablation. However, we knew little about the characteristics of debris cover in this area, except for the proportion of glacier area covered by debris. Moreover, the interactions between glaciers and proglacial lakes, such as ice front calving, ice erosion by lake water, and glacier debris flux into glacial lakes, can also accelerate the glacier retreat and glacial lake expansion. To find out the holistic cause of changes in glacier and glacial lakes, further exploration of above factors needs to be carried out.

Conclusions
In this study, we determined changes in glacial lake area, glacier area and glacier mass balance in the Xixabangma massif in different periods. Our results show that the glacier area shrank from 954.01 km 2 to 752.46 km 2 from 1974 to 2018, whereas the glacial lake area grew from 20.90 km 2 to 38.71 km 2 . Lacustrine-terminating glaciers retreated faster than land-terminating glaciers. During the periods 1974-2000, 2000-2013 and 2013-2018, the region-wide glacier mass balance values were −0.16 m w.e./a, −0.31 m w.e./a and −0.29 m w.e./a, respectively. The annual average temperature during the above three consecutive periods was 3.57 • C, 4.27 • C and 4.22 • C, respectively, indicating a direct response from glacier mass balance to temperature changes. The outburst risk of six glacial lakes was discussed based on the derived changes in lakes and glaciers, as well as other factors. From 1974 to 2018, glacial lakes Gangxico, Galongco and Jialongco expanded by 127.14%, 373.45% and 436.36%, respectively. The mass loss rates of their parent glaciers during 2000-2013 increased by 81.72%, 122.22% and 160.00% relative to losses from 1974 to 2000. We deemed that the outburst risk of glacial lakes Gangxico, Galongco and Jialongco to be high. By contrast, the outburst risk of the other three glacial lakes, Guoqiangco, Guoruco and Kungco, was considered to be low; these latter three lakes pose no threat to local residents and social property.
Priority should be given to detailed field investigations of glacial lakes Gangxico, Galongco and Jialongco, especially regarding the water depth and the distribution of ice cores in lake dams, both of which are difficult to obtain using remote sensing techniques. The terminal condition and water volume of these three glacial lakes should be monitored regularly, and assumptions concerning their risk should be evaluated periodically. Mean-while, it is necessary to install early warning systems for GLOF around high-risk glacial lakes. Enough evacuation time is vital for the purpose of disaster mitigation.