Heterogeneous decadal glacier downwasting at the Mt. Everest (Qomolangma) from 2000 to ~ 2012 based on multi-baseline bistatic SAR interferometry

Remote Sensing of Environment 206 (2018) 336–349

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Heterogeneous decadal glacier downwasting at the Mt. Everest (Qomolangma) from 2000 to ~ 2012 based on multi-baseline bistatic SAR interferometry

T

⁎

Gang Lia, Hui Lina,b,c, , Qinghua Yed a

Institute of Space and Earth Information Science, The Chinese University of Hong Kong, Hong Kong, China Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong, China c Shenzhen Research Institute, The Chinese University of Hong Kong, Shenzhen 518057, China d Institute of Tibetan Plateau Research, Chinese Academic of Science, Beijing 100100, China b

A R T I C L E I N F O

A B S T R A C T

Keywords: Everest Geodetic glacier mass balance Rongbuk TerraSAR-X/TanDEM-X Bistatic D-InSAR

Remote sensing based geodetic observations can be used as alternatives data to map glacier height changes because the harsh environment complicates in situ observations. In this study, we analysed five pairs of X-band bistatic TerraSAR-X/TanDEM-X images from 2011 and 2014 for Mt. Everest (Qomolangma). Glacier height changes were derived using the D-InSAR method respect to SRTM DEM of 2000. An iterative D-InSAR method using multi-baseline bistatic SAR interferograms was employed, which increased reliability and accuracy of glacier height changing observations. From 2000 to ~2012, the geodetic glacier mass balance for the Mt. Everest and the surrounding region was −0.38 ± 0.04 m w.e. a− 1. The spatial pattern of the glacier mass loss was heterogeneous. The regional heterogeneity may possibly reflect debris-covering rates, terminating type, temperature rising rates and glacier flow rates. Comparison to long-period geodetic glacier mass balance data provided by previous studies since 1970 revealed that our results showed more rapid increases in the glacier mass loss rates after 2000 in the area around Khumbu Glacier in the southern slope of the Mt. Everest, whereas glacier mass loss rates kept stable in the Rongbuk Catchment at its northern slope.

1. Introduction The glaciers in the Himalayas and the surrounding areas (High Mountain Asia, HMA) make up the world's largest glacierized area outside the polar regions, and feed the headwaters of many large rivers in Asia (Bolch et al., 2012; Ye et al., 2015). Overall the mass loss from the glaciers in the HMA contributed ~10% of the global glacier mass loss between 2003 and 2009 (Gardner et al., 2013). Affected by the Westerlies, the Indian Monsoon and the East Asian Monsoon, the mass change rates in the different sub-regions of the HMA were heterogeneous (Yao et al., 2012; Gardelle et al., 2013; Gardner et al., 2013; Kääb et al., 2015). Mt. Everest (Qomolangma) is one of the most studied sites in the Central Himalaya, and glacier mass balance observations have been collected for many years (Bolch et al., 2011; Nuimura et al., 2012; Kääb et al. 2012, 2015; Gardelle et al., 2013; Ye et al., 2015). Because of the harsh environment found in mountain environments, field surveys are very challenging tasks to perform, especially at high elevations. The majority of previous studies used satellite laser altimetry or optical

⁎

photogrammetry to calculate glacier height changes (Bolch et al., 2011; Nuimura et al., 2012; Gardelle et al., 2013; Kääb et al. 2012, 2015; Ke et al., 2015; Ye et al., 2015). ICESat footprint data showed that the geodetic glacier mass balance at Mt. Everest and the adjacent region (83°40′E ~ 89°00′E along Himalaya range) during 2003–2008 was − 0.37 ± 0.10 m w.e. (water equivalent) a− 1 (Kääb et al., 2015). Given the large orbital gaps of the ICESat footprints in the HMA, it is impossible to map spatial details at fine scales. Photogrammetry is capable of mapping glacier height changes with better spatial details. A comparison of SPOT5/HRS data obtained in January 2011 with Shuttle Radar Topography Mission (SRTM) data obtained in February 2000 showed that the geodetic glacier mass balance in this region was − 0.26 ± 0.14 m w.e. a− 1 (Gardelle et al., 2013). Another comparison of Worldview 1/2/3 data obtained in 2014 and 2015 to SRTM DEM showed the geodetic mass balance was −0.52 ± 0.14 m w.e. a− 1 (King et al., 2017). Most other studies have examined only a small part of this region with photogrammetry techniques, but spanned longer study periods (Bolch et al., 2011; Nuimura et al., 2012; Ye et al., 2015). Based on a DEM from 1974 and ALOS/PRISM images obtained in 2006,

Corresponding author. E-mail address: [email protected] (H. Lin).

https://doi.org/10.1016/j.rse.2017.12.032 Received 27 October 2016; Received in revised form 29 November 2017; Accepted 27 December 2017 0034-4257/ © 2017 Elsevier Inc. All rights reserved.

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 1. Study site and coverage of the X-band TSX/TDX bistatic SAR images. Polygons represent the TSX/TDX bistatic SAR images coverage. The green polygon shows the TSX/TDX bistatic SAR images pair with a longer perpendicular baseline. The clean ice glacier boundary is identified from Landsat 7 image taken on 30-10-2000. The blue lines divide the study area into sub-regions together with boundaries of Rongbuk catchment (Ye et al., 2015) and Bolch's ROI (Bolch et al., 2011). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

large number of glaciers in the Mt. Everest region using the bistatic differential SAR interferometry (hereafter D-InSAR) method in recent years. To overcome SAR distortion errors, we use cross-track SAR images to fill the foreshortening, layover and shadow regions induced by SAR's side-looking geometry (Li and Lin, 2017). Bistatic image pairs with long and short perpendicular baselines are used to minimize elevation ambiguity (Lachaise et al., 2012a). The relationship between the glacier height and the downwasting rate is analysed, and the geodetic glacier mass balances in different sub-regions are estimated and compared.

the geodetic mass balance of the Rongbuk catchment, on the Chinese side of the Himalayas, was − 0.47 ± 0.23 m w.e. a− 1 (Ye et al., 2015). On the Nepalese side, the average geodetic glacier mass balance of the Khumbu and the surrounding region between 1970 and 2007 was − 0.32 ± 0.08 m w.e. a− 1, and derived rates before and after 2002 indicated an increase in glacier mass loss, though those results were just barely statistically significant (Bolch et al., 2011). Nuimura et al. (2012) studied a larger area on the Nepalese side and obtained a mass loss rate of −0.40 ± 0.25 m w.e. a− 1 from 1992 to 2008. Bistatic SAR interferometry is an alternative remote sensing method for optical photogrammetry and altimetry for analysis of forming topography. It shows potential in forming topography for the glacierized area at Puruogangri ice filed and overcomes the limitation of photogrammetry in observing accumulation zones some extent (Neckel et al., 2013; Liu et al., 2016). The SRTM was launched in February 2000 and was the first satellite mission to use this method to obtain near global topography, with the exception of polar regions (Farr et al., 2007). TanDEM-X (TDX) was launched in 2010 to join its twin satellite TerraSAR-X (TSX) and they operate in bistatic mode. The main goal of this mission was to obtain a global DEM with a vertical accuracy better than 2 m (Gruber et al., 2012). The bistatic observation mode overcomes the temporal decorrelation and atmospheric delay disturbance associated with conventional repeat-pass interferometry (Jaber et al., 2012). However, due to the side-looking geometry, distortion in the foreshortening region can result in some error. In addition, if a region is not fully covered by X-band SRTM data due to its acquisition limitations, then penetration depth difference between the X- and C-band should be evaluated when using C-band SRTM DEM as a reference DEM (Rignot et al., 2001; Hoffmann and Walter, 2006; Gardelle et al. 2012, 2013). This study aims to map the geodetic glacier mass balances for a

2. Study site Mt. Everest (86°55′ E, 27°59′ N) is located in the Central Himalayan chain and along the border between China and Nepal. The main ridge is the national border, which separates Mt. Everest into a northern slope and southern slope, belonging to China and Nepal, respectively. The Rongbuk and Eastern Rongbuk glaciers are located on the northern slope of Mt. Everest within the Rongbuk catchment (Ye et al., 2015). Approximately 6.9% of the glacierized area in Rongbuk catchment is covered by debris (Ye et al., 2015). The debris cover begins to appear at an elevation of ~6100 m on the western branch of the Rongbuk Glacier (also called the Western Rongbuk Glacier), and at ~5550 m on its eastern branch (also called the Mid Rongbuk Glacier). The Khumbu Glacier and its surrounding glaciers are located on the western slope of the ridge between Mt. Everest and Lhotse and have been closely studied in recent years (Bolch et al., 2011; Nuimura et al., 2012; Gardelle et al., 2013; King et al., 2017). Debris cover appears at ~5800 m in the Khumbu region (Nuimura et al., 2012). Randolph glacier inventory (RGI) v5.0 (Arendt et al., 2015) indicates glacier coverage of 337

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

did not change RGI outlines except for several lacustrine-terminating glaciers which presented clear retreating. Several Landsat images were employed to manually identify glacier termini changing (Fig. S1). The clean-ice boundary was identified from the Landsat 7 images taken on 30-10-2000, which presented least could covering rate. To study the geodetic glacier mass balance pattern, we divided glaciers in this region into several sub-regions as shown in Fig. 1 and Table 1. We obtained five acquisitions of bistatic TSX/TDX SAR images in the Co-registered Single Look Slant Range Complex (CoSSC) format from the German Aerospace Center (DLR) in both orbital passes to derive the decadal glacier height changes. These images have resolutions of approximately 2.5 m in both the range and azimuth direction. Long perpendicular baselines can improve the precision of topographic mapping because of their small height of ambiguity. However, this also means large gradient of phase, which makes phase unwrapping difficult. We obtained one pair of images with a long perpendicular baseline (10-10-2012) covering the core area to improve observation precision (green polygon in Fig. 1). The images observed on 04-11-2011 and 0410-2012 were cropped into northern and southern frames when DLR releasing the data. The northern and southern frames were processed separately in interferometric steps and then be mosaiced. Table 2 shows basic information about these bistatic SAR images. C-band SRTM DEM with 30 m resolution (1 arc-second) was chosen as the reference DEM for bistatic D-InSAR processing and as the datum for the unwrapping step since the X-band SRTM DEM does not fully cover the study area due to the narrow swath width. The eastern part of the study site is covered by both C- and X-band SRTM DEM, and this area is used for evaluating penetration depth differences between Cand X-band radar signals (Fig. S1).

Table 1 Information about the glaciers. Bolch's ROI indicates glaciers studied by Bolch et al. (2011). Glacier name

GLIMS ID

Length (m)

Area (Km2)

Sub-region

Unnamed 1 Unnamed 2 Unnamed 3 Kazhenpu Labeilang Lamalangzhuo Unnamed 7 Kangxiong Eastern Group Kada Kadazhanggeli

G087003E27754N G087008E27832N G087020E27939N G087135E27921N G087142E27859N G087111E27962N G087080E27933N G086992E27986N – G087007E28047N G086986E28081N G086985E28103N G086930E28050N G086956E28089N –

8681 13,847 13,879 16,117 11,938 7594 7663 19,178 – 9448 –

17.114 23.845 29.968 46.69 28.275 9.817 13.272 64.393 32.569 19.117 13.093

Southeastern Southeastern Southeastern Southeastern Southeastern Southeastern Southeastern Kangxiong Northeastern Northeastern Northeastern

12,854

39.373

–

80.612

Rongbuk Catchment Northeastern

G086830E28048N

20,118

73.215

G086789E28107N – G086719E28153N –

9882 – 13,185 –

15.916 15.942 47.113 63.412

Rongbuk Catchment Northwestern Northwestern Northwestern Northwestern

G086550E28172N G086557E28122N G086654E28137N G086590E28121N G086536E28082N G086535E28011N –

3677 3422 12,675 17,462 8912 15,952 –

4.079 4.483 19.829 33.291 18.475 27.663 28.354

Northwestern Northwestern Northwestern Western Western Western Western

– G086664E28068N G086716E28089N G086756E28007N G086808E28006N G086903E27975N G086868E27946N G086891E27940N G086951E27927N G086958E27924N G086936E27879N G086885E27873N –

– 10,940 18,285 9530 5056 15,396 5939 4305 5924 8812 2337 4155 –

21.071 12.919 61.054 12.604 12.451 19.097 3.692 2.842 6.825 14.266 1.922 4.824 59.502

Western Western Ngozumpa Ngozumpa Ngozumpa Bolch's ROI Bolch's ROI Bolch's ROI Bolch's ROI Bolch's ROI Southeastern Southeastern Southeastern

Eastern Rongbuk Northeastern Group Rongbuk Jiuda Northern Group Gechongkang Northwestern Group Yanong North Yanong Unnamed 19 Unnamed 20 Shalong Unnamed 22 Southwest Group Western Group Unnamed 25 Ngozumpa Unnamed 27 Changri Khumbu Nuptse Lhotse Nup Lhotse Imja Amphu Laptse Amadablam Southern Group

4. Method 4.1. Bistatic D-InSAR Monostatic interferograms contain both topographic and deformation phase, which can be used to derive glacier height changes, flow velocities, and strain rate by applying two-pass, three-pass, and fourpass methods when the coherence allows (Cheng et al., 2007; Buchroithner and Walther, 2007). One-pass bistatic interferograms in one-pass bistatic contain only the topographic phase (Kubanek et al., 2015). Two methods can be used to derive glacier height changes using bistatic SAR interferometry (Neckel et al., 2013). The first method uses orbital information from bistatic SAR images and reference DEMs to simulate the flat earth and the topographic phase, and it then removes them from the original bistatic interferogram and leaves a differential interferogram; this technique is called the D-InSAR method. The second method involves generating DEMs from bistatic SAR images using standard InSAR methods and then performing common DEM differencing with respect to a reference DEM, which here is the SRTM DEM. In the D-InSAR method, the smaller phase gradients increase the reliability of phase unwrapping since most parts of topographic phase has been simulated and removed while the first method unwraps the total topographic phase (Neckel et al., 2013). The topographic residual phase dominates the differential interferogram in bistatic mode, which is mainly contributed by glacier height changes and can be transformed

1708.2 km2 at our study site (86°18′~87°28′E, 27°36′~28°24′N) and the lowest terminus at 3718 m for Kazhenpu Glacier. Debris covers approximately 27% glacierized region in Central Himalaya, with a higher debris covering rate along the southern slope (Scherler et al., 2011). Supraglacial lakes are common on the debris-covered glaciers and have grown rapidly in both quantity and areas (Ye et al., 2009; Gardelle, et al., 2011). 3. Data Fig. 1 shows the boundary of the RGI v5.0 (Arendt et al., 2015) and the clean-ice boundary. In this region, RGI v5.0 outlines for most glaciers were derived from Landsat images obtained on 30-10-2000. We Table 2 Specification of the bistatic TSX/TDX SAR dataset used. Date

Relative orbit

Orbital pass

Effective perpendicular baseline (m)

Height of ambiguity (m)

Average incidence angle

Master satellite

04-11-2011 31-01-2012 04-10-2012 10-10-2012 11-02-2014

128 128 29 128 29

Ascending Ascending Descending Ascending Descending

113.8 83.2 88.1 195.5 126.2

−49.9 −62.9 71.5 −28.0 53.8

36 34 39 35 41

TSX TSX TSX TSX TSX

338

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

descending orbits and Fig. 3c shows the glacier height changes after producing the image mosaicking.

directly to height changes without forming a new DEM. Therefore, in this study, we mainly employ the D-InSAR method to detect glacier height changes. The basic InSAR and D-InSAR geometry is precisely known (Hanssen, 2001). The principle of using the D-InSAR method to derive height changes from bistatic SAR images with respect to SRTM DEMs is presented in Eq. 1 (Hanssen, 2001). In the bistatic mode, since only one satellite transmits the signal and two satellites receive it, the factor has been changed from 4π to 2π (Krieger et al., 2013; Kubanek et al., 2015).

∆φheight bi = −

2πB⊥ ∆h 2πB⊥ ∆hsrtm 2πB⊥ ∆hresidual ⎞ = −⎛ + λRsinθ λRsinθ λRsinθ ⎠ ⎝

4.2. Estimation of penetration depth difference Penetration depth difference between C- and X-band SAR signals into the glacier surface was estimated by comparing the C- and X-band SRTM (Gardelle et al., 2013). Here we assumed that the similar surface condition in the adjacent area selected for estimating penetration depth differences and our study regions (Fig. S2). Penetration depth difference is supposed to vary with changes in the clean ice/firn/snow density and water content changes, and debris-covered areas also have different penetration depths (Rignot et al., 2001). Since RGI dataset did not provide debris-cover information, we used Landsat images obtained on 30-10-2000 to separate clean ice and debris-covered ice using NDSI (Arendt et al., 2015). Regions with NDSI values > 0.4 that within the RGI boundary were identified as clean ice/firn/snow on glaciers, whereas NDSI values smaller than 0.4 and within the RGI boundary were identified as debris-covered ice (Fig. 4d). Pixels with NDSI value between 0.35 and 0.45 were mostly mixed pixels on the boundaries. We analysed and corrected the penetration depth differences in each 50-m elevation bin. The C-band SRTM DEM with a 30-m resolution was employed in this step because almost no horizontal drift with respect to the X-band SRTM DEM was found. In addition, the void filled region in the C-band SRTM DEM was also masked out by referring to non-void filled SRTM DEM. At the intersections of the ascending and descending orbits of X-band SRTM DEM, datum differences were usually found with respect to adjacent regions that were only covered by an ascending or descending track (Marschalk et al., 2004). We divided the measured area into two parts to mitigate their datum differences (Fig. S2). The penetration depth differences in the off-glacier regions in both parts were assumed to be zero as the acquisition dates of both SRTM and of our obtained TSX/TDX images avoided the main rainy season (Yang et al., 2011; Salerno et al., 2015). The possible snow cover at the offglacier regions should be thin and both C- and X-band should penetrate to the ground (Dehecq et al., 2016).

(1)

where (B⊥) is the perpendicular baseline, (λ) is the wavelength of the radar signal, (θ) is the incidence angle, (R) is geometric distance from the satellite to the scatterer (also called the slant range distance), and (Δh) is height, which can be split into height in SRTM DEM (Δhsrtm) and the height changes (Δhresidual) due to glacier thickening or thinning. The first part of Eq. 1 was simulated with SRTM and orbital information of TSX/TDX images. It assumed that no height change occurred in the offglacier region in RGI v5.0 because the positions of the glacier termini in this region were believed to be stable except lacustrine-terminating glaciers (Bolch et al., 2008; Arendt et al., 2015; King et al., 2017). Because large changes in glacier height (> 50 m) might have occurred during the study period, the use of SRTM DEM to generate simulated SAR images, and lookup tables between the SAR and geographic coordinates might not be very accurate (Quincey et al., 2009). We applied the D-InSAR method iteratively to generate TSX/TDX DEMs and also for overcoming the difficulty in unwrapping interferograms with long perpendicular baselines. In each iteration, for each pair of bistatic SAR images, we took the DEM generated in the last iteration (SRTM DEM for the first iteration) to derive the topographic residual phase. We added the height changes to the reference DEM in each iteration to generate a new DEM for the next iteration, or we considered it as the final TSX/TDX DEM for each frame. We stopped after the third iteration because the topographic residual phase was close to zero. Images acquired on 10-10-2012 cannot be unwrapped correctly due to its long perpendicular baseline in the first iteration, therefore mosaiced DEM was generated from other acquisitions in the first iteration was employed as reference DEM for the second iteration for this acquisition (Lachaise et al., 2012a). Finally, absolute height changes were transferred to annual scale and mosaiced. Fig. 2 presents the wrapped DInSAR interferograms from each iteration for image pairs obtained on 04-10-2012 (the first row, ascending orbit) and 10-10-2012 (the second row, descending orbit). We applied 3 × 3 multi-look processing to the D-InSAR data and 4-times oversampling to the SRTM DEM. A misalignment between two DEMs results in a false ‘hill-shade’ phenomenon in DEM differencing processing (Nuth and Kääb, 2011). To estimate and remove these ‘hill-shades’, we accurately coregistered the simulated and real SAR images in the geocoding step to refine the lookup table between the SAR and geographic coordinates (Neckel et al., 2013). We applied the normalized cross-correlation (NCC) method to image intensity in the frequency domain for coregistration, and an oversampling factor of 2 was applied. No obvious phase differences (in colour) associated with the topographic aspect were identified (Fig. 2c, f & 3), which suggested that the bistatic SAR images were properly aligned to the reference DEM (Nuth and Kääb, 2011). We reestimated the perpendicular baseline and removed the orbital ramp by applying the FFT method to evaluate the D-InSAR interferograms' fringe rates. We employed C-band SRTM v3 DEM without filling in missing data to remove regions in which the SRTM DEM was filled in with other data. Regions in the bistatic TSX/TDX interferograms affected by geometric distortion including foreshortening, layover and shadows cannot be used because of their poor measurement quality. We evaluated and masked out these regions following Li and Lin (2017). Fig. 3a & b show the glacier height changes after masking out geometric distortion area in the region of overlap between the data from the ascending and

4.3. Seasonal effects To evaluate seasonal effects on the glacier height changes, an overlapping region of five bistatic image pairs was selected to compare the TSX/TDX DEM differences. Relative to the data acquired on 201210-10, and after removing the long-term trend, the height differences were determined to be −0.05 m, −0.35 m, 0.71 m, and 0.19 m for the data acquired on 04-11-2011, 31-01-2012, 04-10-2012 and 11-022014, respectively. No obvious trend was found from October to the following February. Meteorological records show that most precipitation occurs in summer on both sides of the Mt. Everest (Yang et al., 2011; Salerno et al., 2015). Annual precipitation at the Dingri meteorological station (20°38′ N, 87°05′ E, 4301 m a.s.l.) on the northern slope was 263 mm a− 1, and was 449 mm in PYRAMID meteorological station (27°58′ N, 86°49′ E, 5035 m a.s.l.) on the southern slope. The precipitation from October to the following April does not exceed 20 mm on the southern slope of Mt. Everest (Salerno et al., 2015). Precipitations on glaciers could be much higher than the off-glacier region. Ice core records in the Eastern Rongbuk glacier suggested an annual accumulation rate of 0.64 m during 1935–2000 (Kaspari et al., 2008). Considering most accumulation and melting occurring in summer while all our topography observation were in the winter season, no seasonal correction was applied. 4.4. Average glacier height changes and converting glacier height changes to geodetic glacier mass balance Topography affects the quality of glacier height change measurements. Due to the steep slopes found at high elevations, glacier height 339

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 2. Wrapped bistatic D-InSAR interferograms from each iteration. (a), (b) and (c) are D-InSAR interferograms based on data acquired on 04-10-2012 for iterations 1, 2 and 3, respectively, and (d), (e) and (f) are D-InSAR interferograms based on data acquired on 10-10-2012 for iterations 1, 2 and 3, respectively.

terminating glaciers using Landsat images obtained on 27-01-2007 (Fig. S1). To convert the glacier height changes to geodetic glacier mass balance, a glacier density of 850 ± 60 kg m− 3 was applied to both the accumulation zone and the ablation zone. The gradient of ELA from north to south makes it difficult to separate the accumulation zone from the ablation zone (Yao et al., 2012; Huss, 2013). In past studies, typically a constant density was applied. Errors in glacier height measurements dominate the final error in geodetic glacier mass balance (Gardelle et al., 2013).

changes were less likely to be measured with acceptable quality in those regions, where voids appeared in the SRTM and TSX/TDX interferograms. Therefore, the average glacier height changes in each 50-m elevation bin were applied to pixels without available measurements. In addition, based on the assumption that elevation changes should be similar at a given elevation, we linearly interpolated our observations of glacier height changes for several elevation bins because they showed random behaviour due to limited measurement points (Gardelle et al., 2013). RGI v5.0 was regarded as glacier outline for most glaciers except for lacustrine-terminating glaciers (Arendt et al., 2015). We applied manually identified outlines for these quick retreating lacustrine-

Fig. 3. Height changes related to SRTM DEM detected by ascending and descending bistatic D-InSAR. (a) Height changes detected with the ascending acquisition (04-11-2011) in the first iteration. The regions affected by foreshortening and layover are mainly on western slopes. (b) Height changes detected with the descending acquisition (04-10-2012) in the first iteration. The regions affected by foreshortening and layover were mainly on western slopes. (c) Average height changes from these two acquisitions in the first iteration.

340

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 4. False colour Landsat images and the clean ice/firn/snow boundary. (a) False colour Landsat 7 images obtained on 30-10-2000. B3 (blue) is shown in red, B4 (NIR) is shown in green, and B5 (SWIR) is shown in red. Clean ice/firn/snow is shown in cyan. (b) Clean ice/firn/snow coverage identified from the Landsat images obtained on 30-10-2000 shows in white. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4.5. Accuracy estimation

Errrandom =

The error in the detected glacier height changes is composed of systematic and random components. The former includes errors in the penetration depth, seasonal effects, and the datum between the SRTM and TSX/TDX data in the off-glacier regions. The random component depends on the precision of the TSX/TDX and SRTM DEMs as well as the total glacier area measured. To estimate the error of the penetration depth, we calculated the standard deviation for each 50-m elevation bin and then take the weighted average as the RMSE of the penetration depth difference. The error in the penetration depth difference was then calculated with Eq. 2: 2

Errpd =

2

σbi − SRTM (On) =

(2)

Errpd 2 + ErrSeason 2 + Errdatum 2

σbi − SRT (Off ) i ∗P(on) i

(6)

5.1. Penetration depth difference and seasonal effect In the Mt. Everest region, the penetration depth difference for clean ice was about 2.26 m below 6000 m, and 2.78 m between 6000 and 6900 m. SRTM-X DEM does not cover Mt. Everest summit and penetration depth difference cannot be determined for regions above 6900 m. According to the linear trend calculated for the several highest bins in Fig. 5 and confirmed by a study on ICESat vs. SRTM (Kääb et al., 2012), we assumed a value of 4.8 m for the area above 6900 m. Glacier area above 6900 m occupies about 2% of the total glacier area in this study. For the debris-covered region, we did not apply any penetration corrections. In total, the average penetration depth difference for C-/Xband SRTM overlapped region was 2.53 ± 0.46 m for clean ice/firn snow.

(3)

where (σbi-SRTM(Off)) is the standard deviation of the height difference in the off-glacier region. (σbi − SRTM) is RMSE of differencing between TSX/TDX data and SRTM, which is composed by the SRTM error, phase noise in the TSX/TDX images and errors induced by slight misalignments between the pairs of DEMs. (Neffbi − SRTM(Off)) is the number of effective measurements, using 500 m as an autocorrelation distance according to the variogram analysis. Then, the systematic component of the error was calculated with Eq. 4:

Errsys =

7

∑i =1

5. Results

σbi − SRTM (Off ) Neffbi − SRTM (Off )

(5)

where (i) is the number of bins. Detailed information about the standard deviation in each slope bin and the proportions of the off-glacier and glacierized areas is listed in Table S1.

where (Errpd) is the error of the penetration depth difference. (σcx_on) is the RMSE of the penetration depth differences (approximately 6.8 m), and (Neffcx_on) is the number of effective measurement on glaciers. (σcx_off) is the RMSE of the penetration depth difference (approximately 6.0 m), and (Neffcx_off) is the number of effective measurement in the off-glacier region. Based on a variogram analysis, the autocorrelation distance was set to 1000 m for the differencing process of the C- and X-band SRTM DEMs. We assumed an error of 0.15 m/month as seasonal component for images not observed in February, which was when the SRTM were obtained (Gardelle et al., 2013). The datum error was then calculated with Eq. 3:

Errdatum =

Neffbi − SRTM (On)

where (Neffbi − SRTM(On)) is the number of effective measurements on glacier, using 500 m as the autocorrelation distance and (σbi-SRTM(On)) is the RMSE of the height difference in the glacier region. However, since this part cannot be measured directly, we assumed that the RMSE of height difference was correlated to the slope rates (Pieczonka et al., 2011). We evaluated the average standard deviations and proportions of glacierized area separately in seven slope bins in the off-glacier region, including 0–5, 5–10, 10–20, 20–30, 30–40, 40–50, and > 50°. Considering the proportion of each bin as weight (P(on)i) and standard deviation (σbi − SRT(Off)i) for each bin in the off-glacier region, the standard deviation for glacierized region was calculated with Eq. 6:

2

⎛ σcx off ⎞ ⎛ σcx on ⎞ ⎟ ⎜ Neff ⎟ + ⎜ Neff cx on ⎠ cx off ⎝ ⎝ ⎠

σbi − SRTM (On)

5.2. Cross validation The inset panels in Fig. 6 present histograms of the height differences in off-glacier regions for each TSX/TDX acquisitions with respect to the SRTM DEM, which was applied as (σbi-SRTM(Off)) in Eqs. 3 and 6.

(4)

The random error is calculated with Formula 5: 341

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 5. Penetration depth differences between C- and X-band SRTM at each elevation bin. Blue indicates clean-ice/firn/snow, while grey indicates debris-covered ice. (a) Penetration depth differences, (b) Glacier height distribution. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

heterogeneous pattern at the sub-regional scale. Fig. 7 shows the glacier height changes and geodetic glacier mass balances. The geodetic glacier mass balance for the Rongbuk catchment (magenta boundary in Fig. 7) between 2000 and ~2012 was − 0.51 ± 0.05 (0.02) m w.e. a− 1, which was more negative than the regional average. The Rongbuk and Eastern Rongbuk glaciers changed more rapidly at −0.69 ± 0.06 (0.02) m w.e. a− 1. The most obvious glacier thinning occurred in the debris-covered areas within the ablation zone where the Mid Rongbuk Glacier joins the Western Rongbuk Glacier at an elevation of ~5480 m. Similar phenomena were observed on several other heavily debris-covered glaciers in this region. The geodetic glacier mass balance of the Khumbu Glacier and its surrounding region (Bolch's ROI; Bolch et al., 2011), was −0.48 ± 0.05 (0.03) m w.e. a− 1 during the period 2000–~2012, whereas the Khumbu Glacier experienced more rapid degradation at − 0.63 ± 0.06 (0.03) m w.e. a− 1. The glaciers in the northeast and northwest sub-regions where clean ice glaciers dominated experienced smaller rates of change than their surrounding glaciers. Most of the glaciers do not extend to much lower elevations, and could be even higher than the middle part of the ablation zone of the Rongbuk and the Eastern Rongbuk glaciers. Several glaciers with large accumulation zones and small ablation zones also have lower glacier mass losses, such as the Gechongkang Glacier and the Kangxiong Glacier. Table 4 provides detailed information about the estimated geodetic glacier mass balances for the total area, the northern and southern sides and each sub-region. The glacier height changes were related to the elevation. Glacier height changes in every 50-m elevation bin for clean ice and debriscovered ice present the results in Fig. 8. For the Chinese side below ~6700 m and the Nepalese side below ~6200 m, the glacier downwasting rates increased nearly linearly with decreasing elevation,

The main panel shows the cross-validation of the glacier height changes detected in the data acquired on 04-11-2011n and 04-10-2012n. The standard deviation of the height difference between these two observations in our study is 3.71 m in the glacierised region. Theoretically, the standard deviation of height differences between pairs of TSX/TDX observations should be 1.414 times the RMSE of the TSX/TDX DEMs, which DLR claimed to be 2 m (Krieger et al., 2013). We also cross-validated the quality of TSX/TDX observations in the overlapped off-glacier regions. RMSEs of height differences are typically better than 3 m, except for several involved with the acquisition on 04-10-2012, which has the largest height ambiguity in the five acquisitions. Detailed information is presented in Table 3. Generally, acquisitions with a smaller height of ambiguity (longer perpendicular baselines) present better precisions than the acquisitions with a larger height of ambiguity (shorter perpendicular baselines).

5.3. Geodetic glacier mass balance Since all the sub-regions have the same error source for the penetration depth difference, we provide the errors of the geodetic glacier mass balance with and without the component of penetration depth difference. The latter should only be used for determining the heterogeneous geodetic glacier mass balance pattern and is shown in the parentheses after the total error. The average geodetic glacier mass balance for Mt. Everest and the surrounding region was − 0.38 ± 0.04 (0.02) m w.e. a− 1. A slight difference in the rates of glacier mass loss was identified between the Chinese and Nepalese sides (northern and southern sides, respectively), which were −0.35 ± 0.04 (0.02) m w.e. a− 1 and − 0.41 ± 0.04 (0.02) m w.e. a− 1, respectively. Although the difference between the northern and southern parts of the Himalaya was not great, the geodetic glacier mass balance still has a 342

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 6. Cross-validation of the glacier height detection. The main panel presents the glacier height changes detected in the data acquired on 04-11-2011 and 04-10-2012. The glacier height changes for the latter were normalized to one decade. The inset panels present histograms of the height changes in off-glacier regions.

Table 3 RMSEs of height difference at the off-glacier regions between each TSX/TDX acquisition. Parentheses after the dates present their height of ambiguity. Date

04-11-2011 (−49.9 m)

31-01-2012 (− 62.9 m)

04-10-2012 (71.5 m)

10-10-2012 (− 28.0 m)

11-02-2014 (53.8 m)

04-11-2011 31-01-2012 04-10-2012 10-10-2012 11-02-2014 Average RMSE

– 2.31 m 3.02 m 1.49 m 2.21 m 2.26 m

2.31 m – 3.92 m 2.53 m 2.52 m 2.82 m

3.02 m 3.92 m – 3.40 m 3.12 m 3.37 m

1.49 m 2.53 m 3.40 m – 2.32 m 2.44 m

2.21 m 2.52 m 3.12 m 2.32 m – 2.54 m

average clean-ice downwasting rates were − 0.49 ± 0.04 m w.e. a− 1 and −0.14 ± 0.04 m w.e. a− 1 for RaER and Bolch's ROI, respectively, whereas the debris-covered parts downwasted at − 1.28 ± 0.05 m w.e. a− 1 and − 0.92 ± 0.05 m w.e. a− 1, respectively. Clearly, the average downwasting rate for the debris-covered ice was higher than that of the clean ice. However, comparison of the melting rates at particular elevations shows that in both RaER and Bolch's ROI glacier thinning rates of the debris-covered ice were generally less than those of clean ice in the same elevation, except for some high elevation bins where the debris started to expose and thin (Fig. 9).

except in very low regions. Above these elevations, the glacier height changes were nearly the same, even with increasing elevation. These levels are much higher than the equilibrium lines, which were between 5800 and 6200 m in the Rongbuk Catchment (Ye et al., 2015) and between 5200 and 5800 m in the Nepalese side (Shea et al., 2015). As these ELAs were recorded a long time ago, ELA could be higher in the present day. The glacier downwasting rates for the clean ice areas were − 0.28 ± 0.04 m w.e. a− 1 and − 0.31 ± 0.04 m w.e. a− 1 on the Chinese and Nepalese sides, respectively, and the debris-covered ice downwasted at − 0.72 ± 0.04 m w.e. a− 1 and − 0.78 ± 0.04 m w.e. a− 1 on the Chinese and Nepalese sides, respectively. Given the clear ELA gradient from north to south (Yao et al., 2012; Ye et al., 2015; Shea et al., 2015), when only one glacier or a small group of glaciers was considered, the effects of debris covering on the glaciers is clear. We take the Rongbuk and Eastern Rongbuk glaciers (RaER) and Bolch's ROI as examples for demonstration (Fig. 9). The

6. Discussion 6.1. Effects of debris cover, ice cliffs and supraglacial lakes Debris cover is typically thought to protect ice from melting when 343

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 7. Glacier height changes and geodetic glacier mass balances. Each value represents the geodetic glacier mass balance in m w.e. a− 1. Errors in parentheses do not contain the error component of penetration depth difference.

and the glacier downwasting. Fig. 10a shows the Rongbuk Lake, and its ice cliff. Fig. 10b shows the position and angle from which Fig. 10a was taken (located by a Nikon GP-1A, and validated by a Holux M241 GPS Logger). Several small supraglacial lakes and their positions from which their photos were taken are also shown in Fig. 10c and 10d. Ice cliffs make important contributions to the total ablation of a debris-covered glacier, although they only occupy a small proportion of the total debris-covered area (Benn et al., 2012; Reid and Brock, 2014). These ice cliffs are too steep to hold large debris, but very fine debris reduced the ice albedo and allowed it to absorb more shortwave radiations, and adjacent warm debris layers also emit longwave radiations (Reid and Brock, 2014). Most of the ice cliffs in the debris-covered part of the Mid-Rongbuk Glacier face north (downstream direction). Fig. 10g is a photo that looks south and shows a large number of ice cliffs where the Mid-Rongbuk Glacier joins the Western Rongbuk Glacier, which is the location of the maximum elevation changes (Fig. 7). For downstream of Rongbuk Lake, although supraglacial lakes were still common, almost no exposed ice cliffs were found, and the slope angles were lower (Fig. 10c & d). The absence of exposed ice cliffs could be the reason that the region downstream of Rongbuk Lake had smaller elevation changes (First two bins in Fig. 9a) than the region upstream of the Rongbuk Lake.

Table 4 Geodetic glacier mass balances and debris cover percentages. The region and sub-region names are as shown in Table 1. The total error and the error without the component of penetration are the values outside and within the parentheses, respectively. Region or sub-region

Geodetic glacier mass balance (m w.e. a− 1)

Debris cover percentage

Total region Chinese (Northern) Nepalese (Southern) Southeastern Kangxiong Northeastern Rongbuk Catchment Northwestern Western Ngozumpa Bolch's ROI

−0.38 −0.35 −0.41 −0.28 −0.13 −0.22 −0.51 −0.19 −0.51 −0.43 −0.48

17.7% 15.3% 20.9% 23.4% 27.0% 3.7% 12.1% 9.0% 31.6% 36.7% 47.2%

± ± ± ± ± ± ± ± ± ± ±

0.04 0.04 0.05 0.05 0.04 0.04 0.05 0.04 0.05 0.05 0.05

(0.02) (0.02) (0.02) (0.03) (0.02) (0.02) (0.02) (0.02) (0.02) (0.02) (0.03)

temperatures increase. Although the debris-covered region had lower glacier thinning rates at most elevations, the sub-regions (Northeastern and Northwestern) with the smallest percentages of debris cover had the lowest glacier loss rates in our study. Other supraglacial features may also affect glacier mass balance. Geodetic mass balance for all lacustrine-terminating glaciers (Rongbuk, Imja, Unnamed 2, Yanong, Yanong North, and Kada) was − 0.51 ± 0.05 (0.02) m w.e. a− 1, which was higher than the regional average (0.38 ± 0.04 (0.02) m w.e. a− 1). Land-terminating glaciers with heavy debris-covers (Labeilang, Kazhenpu, Eastern Rongbuk, Gechongkang, Jiuda, Unnamed 22, 23, 25, Ngozumpa, Khumbu, Changri, and Lhotse) presented average downwasting rate at − 0.39 ± 0.04 (0.02) m w.e. a− 1, which was quite similar to the regional average. Rongbuk Lake increased sharply since 1992 (Ye et al., 2009), and it was almost as wide as the glacier when we visited it in summer of 2015. Supraglacial lakes were common for most land-terminating glaciers, but they were not expanding as quick as the supraglacial lake of the Rongbuk glacier or other proglacial lakes (Ye et al., 2009; King et al., 2017). This suggests a possible correlation between supraglacial/proglacial lake expansion

6.2. Comparison with other results Different glacier density scenarios have been used by different studies, therefore, to compare our results to others, we normalized all the results to an assumed density of 850 ± 60 kg·m− 3 (Huss, 2013). The errors were also normalized to 1σ. Kääb et al. (2015) and Gardner et al. (2013) applied ICESat/GLAS laser altimetry and obtained regional geodetic glacier mass balances of − 0.37 ± 0.10 m w.e. a− 1 and − 0.44 ± 0.10 m w.e. a− 1, respectively for the period of 2003–2009. These two values were similar to the spatial averages near Mt. Everest but were calculated over much larger coverage than in our study. Gardelle et al. (2013) applied SRTM (2000-02) and SPOT/HRS stereo images (2011-01-04). After applying the penetration depth difference between C- and X-band SRTM DEMs at 1.4 m, they concluded that the 344

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 8. Glacier height changes on the Chinese side (a) and the Nepalese side (b). The error bars in this figure indicate the standard deviation of the glacier height change in each elevation bin. Since the glacier height changes in one bin are expected to differ, the standard deviation does not correspond to the standard error in this situation (same for Fig. 9).

average glacier downwasting rate was − 0.26 ± 0.13 m w.e. a− 1. However, if our estimated average penetration depth difference of 2.53 m was applied (see Section 4.2 and Supplementary information), their result was −0.35 ± 0.13 m w.e. a− 1, similar to our results of − 0.38 ± 0.04 m w.e. a− 1. We list the results from previous studies for the different sub-regions in Table 5 and show them in Fig. 11. The results of Gardelle et al. (2013) that were corrected with our estimated penetration depth differences are shown in parentheses. This study and most previous studies did not use the same study regions and periods, which could be a reason for the differences between the estimated geodetic glacier mass balances. Bolch et al. (2011) found increasing glacier mass losses in Khumbu and the adjacent region (Bolch's ROI). They suggested a geodetic glacier mass balances of − 0.31 ± 0.08 m w.e. a− 1 (1970–2007) and −0.75 ± 0.50 m w.e. a− 1 (2002–2007), although the latter observation period was very short (Bolch et al., 2011). Our results gives a higher loss rate (− 0.48 ± 0.05 m w.e. a− 1) than that of Bolch et al. (2011) for the earlier period. Our result for the Rongbuk catchment on the Chinese side (−0.51 ± 0.05 m w.e. a− 1) is almost the same comparing to the result of Ye et al. (2015) for 1974–2006 (−0.47 ± 0.23 m w.e. a− 1).

2011). These results also suggested that the northern slope experienced increasing amounts of precipitation after the 1970s during all four seasons, especially in the spring (Yang et al., 2011), which could explain the less negative mass balances on the northern slope (Chinese side). Glacier downwasting rates were primarily associated with elevations (Figs. 8 & 9). Glacierized area concentrates more at low levels (between 4900 and 5500 m) for Khumbu glacier and surroundings while glacierized in the Rongbuk Catchment are more uniformly distributed in more elevation levels. This also possibly explains the increase in glacier mass loss rate for Khumbu and its surroundings.

6.4. Limitations and potentials Evaluating glacier height changes by comparing DEMs generated by optical photogrammetry and SAR interferometry requires the penetration depths of the SAR signals to be evaluated (Gardelle et al., 2013; Holzer et al., 2015). Usually, it presumes that the X-band signal does not penetrate glaciers, or compare InSAR DEM to laser altimetry observations (Kääb et al., 2012; Gardelle et al., 2013; Holzer et al., 2015; King et al., 2017). The good agreement between the geodetic glacier mass balances obtained by applying the same estimated penetration depth of 2.53 m indicates that it is reasonable to make this assumption in the central Himalaya. However, a recent study found that even Xband can penetrate 4 m into dry snowpack in Mont-Blanc area, which suggest the penetration depth difference should be treated carefully when deriving glacier height changes with DEMs observed by different platforms, especially dry snowpack area occupies great proportion of glacierized area (Dehecq et al., 2016). In our study site, the glacierized area above 6900 m only occupies ~2% of the total glacierized area, therefore large error in estimating penetration depth differences in

6.3. Regional heterogeneity and possible link with climatological variables Meteorological records for the regions between 2660 and 5600 m a.s.l. suggested that the temperature increased at 0.53 °C per decade and that substantial liquid precipitation weakening of − 9.3 mm a− 1 occurred on the southern slope of Mt. Everest between 1994 and 2003 (Salerno et al., 2015). On the northern side, the mean temperature increased by approximately 0.36 °C per decade between 1959 and 2007, which were recorded at the Dingri meteorological station, and more warming occurred in the winter than in the summer (Bolch et al., 345

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 9. Glacier height changes at Rongbuk Glacier, Eastern Rongbuk Glacier (a) and Bolch's group. Bolch's group (Bolch et al., 2011) includes Khumbu and several other glaciers. The error bars indicate the standard deviation of the glacier height changes in each elevation bin.

snowpack area should not contaminate the final result of glacier mass balance. Bistatic SAR interferometry holds great potential for generating DEMs. The standard deviations of the height differences between the SRTM and TSX/TDX DEMs in this study are slightly < 4 m, whereas the comparison of the SRTM and SPOT/HRS DEMs gave a standard deviation of 8.66 m, and comparison of the SRTM and WorldView DEMs gave a standard deviation at 6.11 m (Gardelle et al., 2013; King et al., 2017). In addition, multiple bistatic SAR data sets with different perpendicular baselines could easily identify unwrapping errors because of the different height of ambiguities (Lachaise et al., 2012b). This not only increases the precision in measuring topography but also can spot large errors due to unwrapping. According to our error propagation, the majority of the error in difference processing of SRTM and TSX/TDX DEMs was contributed by the error of the penetration depth difference. The standard error of the estimated penetration depth difference in study region was 0.46 m, and the average standard error of seasonal effects for total study region was 0.20 m. Together with the datum error of 0.07 m, the systematic error was 0.51 m. We evaluated the standard error of the random part of the observations as 0.05 m. Therefore, the standard error of the glacier height change was 0.51 m, which was slightly higher than the standard error of the penetration depth differences; therefore, the error in penetration depth differences estimation dominated the total error in estimating glacier height changes. Bolch's ROI (Bolch et al., 2011) is one of the smallest groups (37 km2) in this study. Our error propagation gave an estimated random part of the error in height change observations of 0.33 m, which was still smaller than the error of penetration depth differences. Therefore, if the penetration depth differences can be

evaluated with higher precision, our estimated annual geodetic glacier mass balances could be more accurate. This suggests that radar signal penetration into glaciers is an active area of research. Topographic differencing between two TSX/TDX DEMs obtained with a certain temporal baseline will likely give a higher measurement accuracy for geodetic glacier mass balance calculations and be capable of providing geodetic glacier mass balance time series because the penetration depth differences can be neglected (Liu et al., 2016). 7. Conclusions In this study, we obtained five pairs of bistatic X-band TerraSAR/ TanDEM-X images with short and long perpendicular baselines observed during 2011 and 2014 that cover Mt. Everest and the surrounding regions and estimated the geodetic glacier mass balances. We proposed an iterative method of applying bistatic images with multi perpendicular baseline to overcome decorrelation and unwrapping errors and to achieve a higher measurement precision for topography. Bistatic images with smaller height of ambiguity (usually longer perpendicular baseline) basically present higher precision in deriving topography. By finding the differences between the C- and X-band SRTM DEMs from 2000, glacier height changes with fine spatial details were mapped in almost all altitude levels. Glacier downwasting rates were primarily associated with elevation. The geodetic glacier mass balance for Mt. Everest and its surroundings from 2000 to ~2012 was − 0.38 ± 0.04 m w.e.a− 1. Glacier degradation was faster on the Nepalese side (− 0.41 ± 0.05 m w.e.a− 1) than on the Chinese side (− 0.35 ± 0.04 m w.e.a− 1), which possibly was linked to the climatological variables. A heterogeneous pattern of glacier downwasting 346

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Fig. 10. Photos for supraglacial features for the Rongkbuk glacier. (a) Rongbuk Lake. (b) Position and approximate view angle of the photo of Rongbuk Lake. (c) Supraglacial lakes. (d) Same as (b) but for (c). (e) Ice cliffs near Rongbuk Lake. (f) Same as (b) but for (e). (g) A photo facing south showing numerous ice cliffs. (h) Same as (b) but for (g). Photos were taken on 31-07-2015.

Acknowledgements

was found in the area surroundings Mt. Everest. Clean ice glaciers and glaciers with high flow rates lost the least glacier mass. Land-terminating debris-covered glaciers' downwasting rate was on par with the regional average. Lacustrine-terminating debris-covered glaciers (including the Rongbuk) present faster mass lost. Compared to two studies that reported geodetic glacier mass balances dating back to the 1970s in the Rongbuk catchment at the Chinese side and the Khumbu Glacier and the surroundings at the Nepalese side, our results show glacier mass loss rates increasing after 2000 for the Khumbu Glacier and its surrounding area, while glacier mass loss rate in the Rongbuk catchment kept stable. The error propagation procedure found that the estimates in the penetration depths differences dominate the final error in estimating geodetic glacier mass balances.

This study is jointly supported by the National Basic Research Program of China (2015CB954103), the Research Grants Council (RGC), General Research Fund of HKSAR, China (CUHK 14233016), National Key R & D Program of China (2017YFA0603103), the National Natural Science Foundation of China (41431070, 41590854, and 41621091). We thank United States Geological Survey (USGS) for providing Landsat optical images and C-band SRTM products freely. Bistatic TerraSAR-X/TanDEM-X SAR images in the CoSSC format were provided by German Aerospace Center (DLR) under project XTI_GLAC6924. X-band SRTM was freely provided by DLR as well. Dr. Adam Devlin refined the manuscript and supplementary information. We also thank four anomalous reviewers who improved this manuscript. 347

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

Table 5 Comparison of the geodetic glacier mass balance results of this study to those from other studies. Please note that study regions having the same name do not exactly overlaps between different studies. Regions

Studies

Methods

Periods

Results

Mt. Everest (Central Himalaya) Mt. Everest (Central Himalaya) Mt. Everest Mt. Everest Mt. Everest Rongbuk Catchment Rongbuk Glacier Rongbuk Catchment Rongbuk Glacier Nepalese Side Nepalese Side Bolch's ROI Bolch's ROI Bolch's ROI Bolch's ROI Ngozumpa Glacier Ngozumpa Glacier

Kääb et al. (2012, 2015) Gardner et al. (2013) Gardelle et al. (2013) King et al. (2017) This study Ye et al. (2015) Gardelle et al. (2013) This study This study Nuimura et al. (2012) This study Bolch et al. (2011) Bolch et al. (2011) Gardelle et al. (2013) This study Gardelle et al. (2013) This study

Laser Altimetry Laser Altimetry SPOT/HRS - SRTM WV - SRTM TSX/TDX - SRTM Photogrammetry SPOT/HRS - SRTM TSX/TDX - SRTM TSX/TDX - SRTM Photogrammetry TSX/TDX - SRTM Photogrammetry Photogrammetry Photogrammetry TSX/TDX - SRTM Photogrammetry TSX/TDX - SRTM

2003–2008 2003–2008 2000–2011 2000–2015 2000–~2012 1974–2006 2000–2011 2000–~2012 2000–~2012 1992–2008 2000–~2012 1970–2007 2002–2007 2000–2011 2000–~2012 2000–2011 2000–~2012

− 0.37 − 0.44 − 0.26 − 0.52 − 0.38 − 0.47 − 0.59 − 0.51 − 0.69 − 0.40 − 0.41 − 0.31 − 0.75 − 0.41 − 0.48 − 0.16 − 0.39

± 0.10 ± 0.10 (− 0.35) ± 0.22 ± 0.04 ± 0.23 (− 0.68) ± 0.05 ± 0.06 ± 0.25 ± 0.05 ± 0.08 ± 0.50 (− 0.50) ± 0.05 (− 0.25) ± 0.05

± 0.13

± 0.16

± 0.21 ± 0.16

Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards. Earth Sci. Rev. 114 (1), 156–174. Bolch, T., Buchroithner, M., Pieczonka, T., Kunert, A., 2008. Planimetric and volumetric glacier changes in the Khumbu Himal, Nepal, since 1962 using Corona, Landsat TM and ASTER data. J. Glaciol. 54, 592–600. Bolch, T., Pieczonka, T., Benn, D.I., 2011. Multi-decadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery. Cryosphere 5, 349–358. Bolch, T., Kulkarni, A., Kääb, A., Huggel, C., Paul, F., Cogley, J., et al., 2012. The state and fate of Himalayan glaciers. Science 336, 310–314. Buchroithner, M.F., Walther, S., 2007. Multiparametric cartographic visualisation of glacier rheology. Cartogr. J. 44 (4), 304–312. Cheng, X., Li, X., Shao, Y., Li, Z., 2007. DINSAR measurement of glacier motion in Antarctic Grove Mountain. Chin. Sci. Bull. 52 (3), 358–366. Dehecq, A., Millan, R., Berthier, E., Gourmelen, N., Trouvé, E., Vionnet, V., 2016. Elevation changes inferred from TanDEM-X data over the Mont-Blanc area: impact of the X-band interferometric bias. IEEE J. Sel. Topics Appl. Earth Observ. in Remote Sens. 9 (8), 3870–3882. Farr, T.G., Rosen, P.A., Caro, E., Crippen, R., Duren, R., Hensley, S., ... Seal, D., 2007. The shuttle radar topography mission. Rev. Geophys. 45 (2). Gardelle, J., Arnaud, Y., Berthier, E., 2011. Contrasted evolution of glacial lakes along the Hindu Kush Himalaya mountain range between 1990 and 2009. Glob. Planet. Chang. 75, 47–55. Gardelle, J., Berthier, E., Arnaud, Y., 2012. Slight mass gain of Karakoram glaciers in the early twenty-first century. Nat. Geosci. 5, 322–325. Gardelle, J., Berthier, E., Arnaud, Y., Kaab, A., 2013. Region-wide glacier mass balances over the Pamir-Karakoram-Himalaya during 1999–2011. Cryosphere 7, 1885–1886. Gardner, A., Moholdt, G., Cogley, J., Wouters, B., Arendt, A., Wahr, J., et al., 2013. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. Science 340, 852–857. Gruber, A., Wessel, B., Huber, M., Roth, A., 2012. Operational TanDEM-X DEM calibration and first validation results. ISPRS J. Photogramm. Remote Sens. 73, 39–49. Hanssen, R.R., 2001. Radar Interferometry. Data Interpretation and Error Analysis. Kluwer Academic Publishers, Dordrecht, The Netherlands. Hoffmann, J., Walter, D., 2006. How complementary are SRTM-X and-C band digital elevation models? Photogramm. Eng. Remote Sens. 72, 261–268. Holzer, N., Vijay, S., Yao, T., Xu, B., Buchroithner, M., Bolch, T., 2015. Four decades of glacier variations at Muztagh Ata (eastern Pamir): a multi-sensor study including Hexagon KH-9 and Pléiades data. Cryosphere 9, 2071–2088. Huss, M., 2013. Density assumptions for converting geodetic glacier volume change to mass change. Cryosphere 7, 877–887. Jaber, W.A., Floricioiu, D., Rott, H., Eineder, M., 2012. Dynamics of fast glaciers in the Patagonia Icefields derived from Terrasar-X and Tandem-X data. In: 2012 IEEE International Geoscience and Remote Sensing Symposium, pp. 3226–3229. Kääb, A., Berthier, E., Nuth, C., Gardelle, J., Arnaud, Y., 2012. Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. Nature 488, 495–498. Kääb, A., Treichler, D., Nuth, C., Berthier, E., 2015. Brief communication: contending estimates of 2003–2008 glacier mass balance over the Pamir–Karakoram–Himalaya. Cryosphere 9, 557–564. Kaspari, S., Hooke, R.L., Mayewski, P.A., Kang, S., Hou, S., Qin, D., 2008. Snow accumulation rate on Qomolangma (Mount Everest), Himalaya: synchroneity with sites across the Tibetan Plateau on 50–100 year timescales. J. Glaciol. 54 (185), 343–352. Ke, L., Ding, X., Song, C., 2015. Heterogeneous changes of glaciers over the western Kunlun Mountains based on ICESat and Landsat-8 derived glacier inventory. Remote Sens. Environ. 168, 13–23. King, O., Quincey, D.J., Carrivick, J.L., Rowan, A.V., 2017. Spatial variability in mass loss of glaciers in the Everest region, central Himalayas, between 2000 and 2015. Cryosphere 11 (1), 407–426. Krieger, G., Zink, M., Bachmann, M., Bräutigam, B., Schulze, D., Martone, M., 2013.

Fig. 11. Comparison of the geodetic glacier mass balances for the entire region and subregions around Mt. Everest. Different icons represent different regions. Colours of the error bars represent different studies. Black: this study, all horizontal; blue: Kääb et al. (2015); green: Gardner et al. (2013); red: Gardelle et al. (2013); magenta: Ye et al. (2015); yellow: Nuimura et al. (2012); cyan: Bolch et al. (2011); orange: King et al. (2017). For Khumbu glacier, Bolch et al. (2011) provided a time series of geodetic glacier mass balance, we adopt the period of 1970–2007 with the minimum error here to compare with our result. Geodetic glacier mass balance values for Gardelle et al.'s (2013) research were corrected estimated penetration differences from their 1.4 m to our estimation of 2.53 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Appendix A. Supplementary data Supplementary data to this article can be found online at https:// doi.org/10.1016/j.rse.2017.12.032. References Arendt, A., Bliss, A., Bolch, T., Cogley, J.G., Gardner, A.S., Hagen, J.-O., et al., 2015. Randolph glacier inventory – a dataset of global glacier outlines: version 5.0. Digital media, Boulder Colorado, USA. In: Global Land Ice Measurements from Space. Benn, D.I., Bolch, T., Hands, K., Gulley, J., Luckman, A., Nicholson, L.I., et al., 2012.

348

Remote Sensing of Environment 206 (2018) 336–349

G. Li et al.

ISPRS J. Photogramm. Remote Sens. 66, 927–940. Quincey, D.J., Luckman, A., Benn, D., 2009. Quantification of Everest region glacier velocities between 1992 and 2002, using satellite radar interferometry and feature tracking. J. Glaciol. 55 (192), 596–606. Reid, T.D., Brock, B.W., 2014. Assessing ice-cliff backwasting and its contribution to total ablation of debris-covered Miage glacier, Mont Blanc massif, Italy. J. Glaciol. 60 (219), 3–13. Rignot, E., Echelmeyer, K., Krabill, W., 2001. Penetration depth of interferometric synthetic-aperture radar signals in snow and ice. Geophys. Res. Lett. 28, 3501–3504. Salerno, F., Guyennon, N., Thakuri, S., Viviano, G., Romano, E., Vuillermoz, E., et al., 2015. Weak precipitation, warm winters and springs impact glaciers of south slopes of Mt. Everest (central Himalaya) in the last 2 decades (1994–2013). Cryosphere 9, 1229–1247. Scherler, D., Bookhagen, B., Strecker, M.R., 2011. Spatially variable response of Himalayan glaciers to climate change affected by debris cover. Nat. Geosci. 2011 (4), 156–159. Shea, J.M., Immerzeel, W.W., Wagnon, P., Vincent, C., Bajracharya, S., 2015. Modelling glacier change in the Everest region, Nepal Himalaya. Cryosphere 9 (3), 1105–1128. Yang, X., Zhang, T., Qin, D., Kang, S., Qin, X., 2011. Characteristics and changes in air temperature and glacier's response on the north slope of Mt. Qomolangma (Mt. Everest). Arct. Antarct. Alp. Res. 43, 147–160. Yao, T., Thompson, L., Yang, W., Yu, W., Gao, Y., Guo, X., et al., 2012. Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings. Nat. Clim. Chang. 2, 663–667. Ye, Q., Zhong, Z., Kang, S., Stein, A., Wei, Q., Liu, J., 2009. Monitoring glacier and supraglacier lakes from space in Mt. Qomolangma region of the Himalayas on the Tibetan Plateau in China. J. Mt. Sci. 6, 211–220. Ye, Q., Bolch, T., Naruse, R., Wang, Y., Zong, J., Wang, Z., et al., 2015. Glacier mass changes in Rongbuk catchment on Mt. Qomolangma from 1974 to 2006 based on topographic maps and ALOS PRISM data. J. Hydrol. 530, 273–280.

TanDEM-X: a radar interferometer with two formation-flying satellites. Acta Astronautica 89, 83–98. Kubanek, J., Westerhaus, M., Schenk, A., Aisyah, N., Brotopuspito, K.S., Heck, B., 2015. Volumetric change quantification of the 2010 Merapi eruption using TanDEM-X InSAR. 164, 16–25. Lachaise, M., Balss, U., Fritz, T., Breit, H., 2012a. The dual-baseline interferometric processing chain for the TanDEM-X mission. In: 2012 IEEE International Geoscience and Remote Sensing Symposium, pp. 562–5565. Lachaise, M., Fritz, T., Balss, U., Bamler, R., Eineder, M., 2012b. Phase unwrapping correction with dual-baseline data for the TanDEM-X mission. In: 2012 IEEE International Geoscience and Remote Sensing Symposium, pp. 5566–5569. Li, G., Lin, H., 2017. Recent decadal glacier mass balances over the Western Nyainqentanglha Mountains and the increase in their melting contribution to Nam Co Lake measured by differential bistatic SAR interferometry. Glob. Planet. Chang. 149, 177–190. Liu, L., Jiang, L., Sun, Y., Yi, C., Wang, H., Hsu, H., 2016. Glacier elevation changes (2012–2016) of the Puruogangri ice field on the Tibetan Plateau derived from bitemporal TanDEM-X InSAR data. Int. J. Remote Sens. 37, 5687–5707. Marschalk, U., Roth, A., Eineder, M., Suchandt, S., 2004. Comparison of DEMs derived from SRTM/X-and C-Band. In: 2014 IEEE International Geoscience and Remote Sensing Symposium, pp. 4531–4534. Neckel, N., Braun, A., Kropáček, J., Hochschild, V., 2013. Recent mass balance of the Purogangri Ice Cap, central Tibetan Plateau, by means of differential X-band SAR interferometry. Cryosphere 7, 1623–1633. Nuimura, T., Fujita, K., Yamaguchi, S., Sharma, R., R., 2012. Elevation changes of glaciers revealed by multitemporal digital elevation models calibrated by GPS survey in the Khumbu region, Nepal Himalaya, 1992–2008. J. Glaciol. 58, 648–656. Nuth, C., Kääb, A., 2011. Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change. Cryosphere 5, 271–290. Pieczonka, T., Bolch, T., Buchroithner, M., 2011. Generation and evaluation of multitemporal digital terrain models of the Mt. Everest area from different optical sensors.

349