Next Article in Journal
From Space to the Rocky Intertidal: Using NASA MODIS Sea Surface Temperature and NOAA Water Temperature to Predict Intertidal Logger Temperature
Next Article in Special Issue
Preface: Land Surface Processes and Interactions—From HCMM to Sentinel Missions and Beyond
Previous Article in Journal
Seasonal Change in Wetland Coherence as an Aid to Wetland Monitoring
Previous Article in Special Issue
Evaluation of Methods for Aerodynamic Roughness Length Retrieval from Very High-Resolution Imaging LIDAR Observations over the Heihe Basin in China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assessing Orographic Variability in Glacial Thickness Changes at the Tibetan Plateau Using ICESat Laser Altimetry

1
Department of Geomatics Engineering, Ho Chi Minh City University of Technology, 703500 Ho Chi Minh City, Vietnam
2
Department of Geoscience and Remote Sensing, Delft University of Technology, 2628 CN Delft, The Netherlands
*
Author to whom correspondence should be addressed.
Submission received: 29 September 2016 / Accepted: 9 February 2017 / Published: 15 February 2017

Abstract

:
Monitoring glacier changes is essential for estimating the water mass balance of the Tibetan Plateau. In this study, we exploit ICESat laser altimetry data in combination with the SRTM DEM and the GLIMS glacier mask to estimate trends in change in glacial thickness between 2003 and 2009 on the whole Tibetan Plateau. Considering acquisition conditions of ICESat measurements and terrain surface characteristics, annual glacier elevation trends were estimated for 15 different settings with respect to terrain slope and roughness. In the end, we only included ICESat elevations acquired over terrain with a slope below 20° and a roughness at the footprint scale below 15 m. With this setting, 90 glaciated areas could be distinguished. The results show that most of observed glaciated areas on the Tibetan Plateau are thinning, except for some glaciers in the northwest. In general, glacier elevations on the whole Tibetan Plateau decreased at an average rate of −0.17 ± 0.47 m per year (m a−1) between 2003 and 2009, taking together glaciers of any size, distribution, and location of the observed glaciated area. Both rate and rate error estimates are obtained by accumulating results from individual regions using least squares techniques. Our results notably show that trends in glacier thickness change indeed strongly depend on the relative position in a mountain range.

Graphical Abstract

1. Introduction

The Tibetan Plateau has steep and rough terrain and contains ~37,000 glaciers, occupying an area of ~56,560 km2 [1]. Recent studies report that the glaciers have been retreating significantly in the last decades. The magnitude of glacial change in the last 30 years is location dependent, with the largest reduction in glacial length and area occurring in the Himalayas (excluding the Karakoram) [2]. In the Tien Shan Mountains in the northwest of the Tibetan Plateau, glacier shrinkage has also occurred during the period between 1950 and 2000 [3]. In the Qilian Mountain Region, 910 glaciers have rapidly reduced in area between 1956 and 2003, with a mean reduction of 0.10 km2 per individual glacier, corresponding to a mean rate of 2127 m2 a−1 [4], or a shrinkage of the total glacier area by 30% ± 8% between 1956 and 2010 [5]. In the western Nyaiqentanglha Range, the glacier area decreased by −6.1% ± 3% between 1976 and 2001 [6]. Additionally, the total glacier areas in the inner Tibetan Plateau and in the Himalayas have also retreated between 1970s and 2000s [7,8,9]. Most of the above results were estimated using topographic maps, in situ measurements, and optical remotely sensed images. In recent years, remote sensing techniques such as photogrammetry, interferometry and radar and laser satellite altimetry have been used for assessing vertical glacial and ice-sheet change both on and off the Tibetan Plateau.
Regional changes in ice elevation in the central Karakoram were obtained by determining the difference between two Digital Elevation Models (DEMs), one obtained from the 2000 Shuttle Radar Topographic Mission (SRTM), and one constructed from Satellite Pour l’Observation de la Terre (SPOT5) optical images obtained in 2008 [10]. Advantage of using full DEMs like in this study is that the complete area of interest is covered and that change results are interpretable down to the resolution of the used DEMs. Availability of such DEMS is limited however, therefore it is almost impossible to estimate trends in elevation variation from full DEMs only [11].
Additional elevation data of high quality but sparse spatial coverage were obtained by the Ice, Cloud and land Elevation Satellite (ICESat) mission between 2003 and 2009. Primary mission goal was to study ice sheet mass balance over polar areas [12], but in recent years also ICESat Geoscience Laser Altimeter System (GLAS) data have been exploited to monitor glaciers in mountain regions such as Himalayas, Alps and the Tibetan Plateau. Using ICESat data for change assessment has several challenges. ICESat only sampled elevations along track, while adjacent tracks are separated by ~70 km on the rough Tibetan Plateau; ICESat tracks are repeated but only in an approximate sense: in general two consecutive tracks of the same orbit have non-overlapping ground footprints; moreover, ICESat was not measuring continuously, but only in 18 campaigns of ~1 month resulting in data of different quality due to variations in the laser power.
In [13], ICESat measurements were combined with the SRTM 2000 DEM to obtain glacial thinning and thickening trends over the Himalayas. The SRTM DEM data were not only used as direct reference elevations but also contributed in the design of criteria on, e.g., slope in combination with Landsat data that were used to assess location and state of the glaciers. The sparse sampling of the ICESat data required to regionally group and average available ICESat elevations resulting in regional change trends. In [14], a similar approach of comparing ICESat elevations to a reference DEM was compared to direct differencing between almost repeated along-track ICESat elevations over the European Alps. A worldwide analysis of glacial change based on ICESat elevations, partly extending the results in [13], and compared to mass change estimates from Gravity Recovery and Climate Experiment (GRACE) satellite gravimetry measurements, is given in [15].
Glacial thickening and thinning estimates on the Tibetan Plateau based on ICESat measurements were first reported on in [16]. Again, the SRTM 2000 DEM was used as a basis to obtain ICESat elevation changes. The results indicated that most of the glacial sub-regions had a negative trend in glacial thickness change, excluding one sub-region in the Western Kunlun Mountains in the northwest of the Tibetan Plateau. The sampled glacial sub-regions in [16] were relative large. Consequently, we consider their glacial conditions as not being homogeneous, due to, e.g., orographic precipitation and variation in solar radiation. The significant influence of climatic parameters and spatial variability on glacial change rates has already been demonstrated for several individual glaciers on the Tibetan Plateau [6,17]. In addition, the quality of ICESat elevations is known to be strongly dependent on terrain characteristics.
Therefore, in this study we exploit, as in [16], ICESat/GLAS data for monitoring glacial thickness changes on the whole Tibetan Plateau. Compared to [16] however, our division in glaciated areas is completely different, as in our study we notably incorporate the glacier orientation in the design of regions. In particular, in our study, mountain ranges are divided into two contrasting regions as we separate the ranges with respect to the main center ridge. This allows us to demonstrate the expected effects of orographic precipitation and variations in solar radiation. To do so, an explicit comparison of our results to those in [16] is included in the Discussion Section. In addition, we explore the ICESat/GLAS data by applying different criteria impacting the quality of footprints including acquisition condition and terrain surface characteristics.

2. Materials and Methods

In this section, we describe input elevation data and glacier outlines. Then, we define and build a dataset for monitoring glacier elevation changes. Finally, we clean the dataset and estimate temporal elevation trends of sampled glaciers on the Tibetan Plateau.

2.1. Input Data

Main data sources used to estimate glacier elevation changes at the Tibetan Plateau consist of ICESat/GLAS data, the Global Land Ice Measurements from Space (GLIMS) glacier mask and the Shuttle Radar Topography Mission digital elevation model (SRTM DEM). The ICESat/GLA14 data support land surface elevations between 2003 and 2009. The GLIMS glacier outlines represent the glacial regions on the Tibetan Plateau. The SRTM data show land surface elevations in 2000, used as a base map to be compared with later elevations derived from the ICESat/GLA14 data. To integrate them, all these data are projected onto the World Geodetic System 1984 (WGS84) in horizontal and the Earth Gravitational Model 2008 (EGM2008) in vertical.

2.1.1. ICESat/GLA14 Data

The ICESat/GLAS products are provided by the National Snow and Ice Data Center (NSIDC). Here, we exploit the level-2 GLA14 data, supporting global land surface altimetry between 2003 and 2009 [18]. The GLA14 data are distributed in binary format and are converted into ASCII columns by the NSIDC GLAS Altimetry elevation extractor Tool (NGAT). The geospatial accuracy of each footprint is reported as 5 m in horizontal and 10 cm in vertical for slopes below 1° [19]. The vertical accuracy is strongly dependent on terrain characteristics. In this study, necessary measurements of each footprint extracted from the GLA14 data consist of acquisition time, latitude, longitude, elevation above WGS84, EGM2008 geoid height, saturation correction flag, and number of peaks. The saturation correction flag indicates if elevation data were possibly affected by saturation effects. The number of peaks in the Gaussian waveform decomposition directly relates to land surface geometry [20]. For each ICESat campaign, the ASCII data are converted into the GIS shapefile format, using the location of each footprint. Figure 1 shows the ICESat L2D campaign tracks from 25 November to 17 December 2008 crossing over the Tibetan Plateau.

2.1.2. SRTM DEM

The Shuttle Radar Topography Mission was flown in February 2000 and collected the first ever high resolution near-global digital elevation data. In this study, we use the SRTM 90 m DEM, produced by NASA [21]. This DEM has a resolution of 90 m at the equator corresponding to 3-arc seconds and is distributed in 5 × 5 tiles. To cover the full Tibetan Plateau, 20 SRTM DEM tiles are concatenated, as shown in Figure 1. The tiles are available in both ArcInfo ASCII and GeoTiff format. The digital elevation data were stored in a grid as m × n matrix. The data are projected in a Geographic (latitude/longitude) projection, with the WGS84 horizontal datum and the EGM96 vertical datum. The vertical error of the DEMs is reported to be less than 5 m on relative flat areas and 16 m on steep and rough areas [22]. Note that meanwhile the United Stated Geological Survey (USGS) released a 1 arc second version of SRTM [23]. This version could not yet be incorporated in this study however. In a study over Swiss alpine glaciers [24], SRTM 90 m DEM elevations were compared to the Swiss 25 m national DEM (DHM25). Results indicate that some significant differences occur over individual glaciers, but that differences tend to level out when averaged over larger areas. As we also upscale results from individual glaciers to regions, we may assume that our results obtained using the 90 m SRTM DEM are still valid. Future studies are recommended to check and if necessary correct DEM data for possible co-registration errors [25].

2.1.3. GLIMS Glacier Outlines

The GLIMS project is a project designed to monitor the world’s glaciers, primarily using data from optical satellite instruments. Now, over 60 institutions worldwide are involved in GLIMS for inventorying the majority of the world’s estimated 160,000 glaciers. These glaciers are distributed in GIS shapefile format and are referenced to the WGS84 datum. In this study, we downloaded the glacier mask presenting glacial outlines on the Tibetan Plateau, submitted by Chinese Academy of Sciences, as shown in Figure 1 [1]. The glacier mask is based on aerial photography, topographic maps and in situ measurements. The product was released on 21 July 2005, but the state of the glaciers is expected to represent the situation in 2002 [26]. Each glacier is represented by a polygonal vector with attributes such as identification code, area, width, length, min elevation, max elevation, and name. Note that a new Chinese glacier inventory was published in 2015 [27]; this version could not yet be incorporated in this study.

2.2. Methods

To estimate a glacial thickness change trend, we consider differences between glacial surface elevations derived from 2003–2009 ICESat laser altimetry and a digital elevation model. Here, the digital elevation model is used as a reference surface. In addition, a glacier mask is used to identify ICESat elevations that are likely to sample glaciers. Each difference is time-stamped by the ICESat acquisition time. Valid differences obtained during the same ICESat campaign track over a certain homogeneous glaciated area, also called a sampled glaciated area, are used to estimate a mean difference. Mean differences for each sampled glaciated area are grouped to form a time series. Consecutively, a temporal trend is estimated through the mean differences per area, resulting in a temporal trend of glacial thickening or thinning.
Additionally, differences between the ICESat GLA14 elevations and the reference SRTM DEM may correspond to change in glacial thickness between 2003 and 2009 if certain requirements are met. However, the vertical accuracy of each ICESat footprint strongly depends on terrain surface characteristics, so we have to remove uncertain footprints before the estimation. Therefore, we estimate surface slope and roughness from the SRTM DEM.

2.2.1. Estimating Surface Slope and Roughness from SRTM DEM

Based on the SRTM DEM, the terrain surface parameters slope S and roughness R are estimated, using a 3 × 3 kernel scanning over all pixels of the grid. For each pixel, the slope S in decimal degrees is locally estimated by Equations (1)–(3) [28,29].
S = 180 π × a r c t a n ( d z d x ) 2 + ( d z d y ) 2
d z d x = ( h 3 + 2 × h 6 + h 9 ) ( h 1 + 2 × h 4 + h 7 ) 8 × Δ l o n
d z d y = ( h 7 + 2 × h 8 + h 9 ) ( h 1 + 2 × h 2 + h 3 ) 8 × Δ l a t
where the hi values (i = 1 ÷ 9) are corresponding to the DEM elevations in the kernal while Δlat and Δlon are the width and the height of a grid cell in meters, estimated by distance Equation (4) [30].
d = r × 2 × a tan 2 ( a , 1 a ) a = sin 2 ( φ 2 φ 1 2 ) + cos ( φ 1 ) × cos ( φ 2 ) × sin 2 ( λ 2 λ 1 2 )
where d is the shortest distance over the earth’s surface—the “as-the-crow-flies” distance—between the two points (λ1, ϕ1) and (λ2, ϕ2) in radians in a geographic coordinate system and r is the earth’s radius (mean radius = 6371 km).
The roughness R in meters is defined as the root mean square of the differences e ^ i between the grid heights and the local 3 × 3 plane, best fitting in the least squares sense [31], following Equation (5).
R = i = 1 i = 9 e ^ i 2 9

2.2.2. Determining a Sampled Glaciated Area

Because of the orbital configuration of ICESat and its along track only sampling, Tibetan glaciated areas are only sampled sparsely by ICESat. In addition, surface elevation changes on these mountain glaciers are expected to be affected significantly by the orientation and face of the corresponding mountain range. For example, the south face of the Himalayas is experiencing more precipitation than the north face, while on the other hand north faces experience less incoming solar radiation. Therefore we decided to group nearby glaciers having similar orientation into one sampled glaciated area while, on the other hand, glaciers on different sides of a mountain range ridge were grouped into different areas. First, we extracted footprints of all ICESat campaigns within the GLIMS glacier outlines, as illustrated in Figure 2. Then, each glaciated area outline was manually determined, by considering the locations of the glaciers and the ICESat footprints. For example, in Figure 2 the ICESat-sampled glaciers having a northern orientation were grouped into one glaciated area, A, while those on the other side of the mountain ridge were grouped into another glaciated area, B. Finally each glaciated area was coded by an identification number.

2.2.3. Identifying a Glacier Elevation Difference

A glacier elevation difference Δh is identified as the difference between an elevation of an ICESat footprint within a sampled glaciated area and the reference SRTM DEM, compare Equation (6), where Δh is in meters above EGM2008.
Δh = hICESat − hSRTM = (Elev − GdHt) − (SRTM_elev + 96_08_Ht)
Each glacier elevation difference Δh depends on the characteristics of the terrain illuminated by the ICESat pulse and the characteristics of the ICESat measurement itself. It is in principle also affected by the local quality of the SRTM reference elevation, but in this study it is assumed that the quality of the STRM DEM is not location dependent. What is assessed in this study is the quality of the elevation difference with respect to the attributes described in Table 1. For this purpose, we extract ICESat footprints within the sampled glaciated areas and obtain their full attributes.
A glacier elevation difference Δh is maintained for further analysis if the corresponding ICESat measurement is considered good according to the following criteria. First we select those footprints whose return echo is not or only lightly saturated and moreover have only one peak in its Gauss decomposition. That is the value of SatFlg should equal 0 or 1, and the value of NumPk should equal 1. A footprint with one mode is expected to correspond to homogeneous land surface. Then we remove footprints affected by clouds. If ICESat footprints are affected by clouds, the elevation variation within one track can be very large, while the altitude difference with other tracks is high [16]. In this study, if the ICESat elevation difference to the SRTM DEM Δh is larger than 100 m, the footprint is assumed to be affected by clouds and removed from further analysis.

2.2.4. Different Settings with Respect to Slope and Roughness

Here, we analyze different settings incorporating the terrain surface characteristics slope and roughness. We remove footprints with a slope S bigger than a threshold S0 and roughness R bigger than a threshold R0. Applying strict thresholds will result in a relative small number of remaining glacier elevation differences albeit of relatively high quality. A slope S below 10° is always considered good while a slope of over 30° results in an inacceptable bias. The roughness R is estimated directly from the SRTM data, its lower limit of 5 m corresponds to relative flat areas while its upper limit of 15 m corresponds to high relief and rough areas. In the following, we consider 15 different settings with slope and roughness values within these outer limits, as described in Table 2. Each record in Table 2, corresponding to one such setting, also summarizes the corresponding resulting trend in glacial thinning/thickening for the whole Tibetan Plateau between 2003 and 2009, as determined by the following steps.

2.2.5. Obtaining Mean Glacier Elevation Differences

For each sampled glaciated area, glacier elevation differences all are time-stamped by ICESat acquisition time. The ICESat acquisition time ti is defined per ICESat track segment, where one track is sampling a glaciated area with consecutive individual footprints. A mean glacier elevation difference Δ h ¯ i is considered representative for the height of the glaciated area above the SRTM base map at ICESat acquisition time ti. The mean difference Δ h ¯ i and its standard deviation si is computed using Equations (7) and (8), where k is the number of ICESat footprints in the track segment that are sampling a glaciated area at ICESat acquisition time ti and Δ h i j is the jth elevation difference, j = 1 ÷ k.
Δ h i ¯ = 1 k j = 1 j = k Δ h i j
s i = 1 k j = 1 j = k ( Δ h i j Δ h i ¯ ) 2
Each ICESat acquisition time ti is considered as an epoch in the time series used to estimate a temporal trend using linear regression. Here we only use the mean glacier elevation difference Δ h ¯ i in a time series if its standard deviation si is less than a threshold Std0 and the number of ICESat footprints k is at least six footprints. The threshold Std0 is defined to be equal to the roughness threshold R0 for each setting with respect to terrain slope and roughness. To remove unreliable elevation differences, we build an iterative algorithm. That is, if si is bigger than Std0 and | Δ h i j Δ h ¯ i | is maximal for j in 1 ÷ k, the jth elevation difference Δ h i j is removed. Then, Δ h ¯ i and si are re-computed. This process is repeated until si drops below Std0 or k is less than six. In Figure 3, the values Δ h ¯ i and si representing mean glacier elevation differences and their standard deviations are shown between 2003 and 2009 for two glaciated areas A and B in case that S0, R0, and Std0 are 15°, 10 m, and 10 m, respectively.

2.2.6. Estimating a Temporal Glacial Thickness Change Trend

For each glaciated area on the Tibetan Plateau, a temporal linear trend is estimated if there are at least six average differences or epochs available, corresponding to at least six ICESat campaign tracks during the observed period 2003–2009. For example, Figure 3 shows the distribution of the average differences of the glaciated areas A and B between 2003 and 2009. An annual glacial thickness change trend is estimated by linear adjustment using Equation (9) [33].
x ^ = ( A T A ) 1 A T y
where,
  • y = [ Δ h 1 ¯ Δ h 2 ¯ Δ h n ¯ ] T : the vector of the average elevation differences per epoch.
  • x ^ = [ x 0 v ] : the vector of parameters of the linear trend, offset x0 and rate v.
  • A = [ 1 1   1 t 1 t 2   t n ] T : the design matrix, in which ti denotes the ith epoch.
Note that n is required to be at least six epochs. The rate v of a linear glacial thickness change is obtained by solving Equation (9) and the root mean square error (RMSE), as standard deviation of residuals, is also computed, using Equation (10) with the least-square residual vector e ^ = y A x ^ . This value consists of a combination of possible data errors and mainly the non-validity of the linear regression model.
R M S E = i = 1 i = n e ^ i 2 n
In addition, the propagated standard deviation σvv of the estimated rate v is given in Equation (11), where Q y y denotes the variance matrix, in which si is the standard deviation of the ith average difference. These values are considered as the confidence interval for the estimated glacial thickness change.
Q x ^ x ^ = [ σ x 0 x 0 2 σ x 0 v 2 σ v x 0 2 σ v v 2 ] = ( A T Q y y 1 A ) 1
Q y y = [ s 1 2 0 0 s 2 2 0 0 0 0 s n 2 ]
Continuing to the example of Figure 3, glaciated area A has an elevation decrease of −1.66 ± 0.42 m a−1 and a RMSE of 3.46 m while glaciated area B has an elevation increase of 0.50 ± 0.31 m a−1 and a RMSE of 3.37 m between 2003 and 2009.

3. Results

Following the method above, temporal glacial thickness change trends on the whole Tibetan Plateau between 2003 and 2009 are estimated for 15 different settings with respect to terrain slope and roughness. The results are shown in Table 2. It indicates that, as expected, the number of observed glaciated areas and the RMSEs of differences estimated by the linear regression increase if the thresholds on slope S0 and roughness R0 are relaxed. In practice, the mean rates of glacial thickness change trends on the whole Tibetan Plateau for the five settings from S11 to S15 (all with R0 = 15 m) are quite similar. In addition, the number of trends having a RMSE of over 5 m significantly increases when ICESat footprints at slopes of over 20° are incorporated as well. A RMSE of over 5 m could correspond to a large fluctuation in glacial thickness or a bad fit of the linear trend model.
Here, S0 and R0 are terrain slope and roughness thresholds, respectively. For each setting, N is the number of glaciated areas observable with a given setting and the numbers v ¯ and σ ¯ v v are the resulting overall rate and its propagated standard deviation of glacial thickness change while R M S E ¯ is the average of the root mean square errors (RMSEs) of the linear regression model. Additionally, N5 is the number of observed glaciated areas having a RMSE of below 5 m.
In this paper, we present the results of setting S13, where S0 and R0 equal 20° and 15 m, respectively, because in this case, a maximum number of 67 areas are observed with RMSE ≤ 5 m. We assume that ICESat footprints selected for estimation of glacial thickness change given these settings are relatively appropriate given the steep and rough terrain of the Tibetan Plateau and given the quality of the SRTM DEM.

3.1. Overall Glacial Thickness Changes: Tibetan Plateau and Its Basins

In the case the thresholds S0 = 20° for terrain slope and R0 = 15 m for roughness are applied, the result indicates that 90 glaciated areas on the whole Tibetan Plateau are sampled by enough ICESat footprints to estimate thickness change. In addition, 67 RMSEs are below 5 m. For each glaciated area, a temporal trend in glacial thickness is estimated, as shown in Table S1. In Figure 4, a glacial thickness change rate is symbolized by a red or blue disk at a representative location in each observed glaciated area. Most of the observed glaciated areas in the Himalaya, the Hengduan Mountains and the Tanggula Mountains experienced a serious decrease in glacial thickness. However, in most of the observed glaciated areas in the Western Kunlun Mountains in the northwest of the Tibetan Plateau, glaciers oriented toward the north were thickening while those oriented toward the south were thinning. In general, glacial thickness on the whole Tibetan Plateau decreased between 2003 and 2009 at a mean rate of −0.17 ± 0.47 m a−1. This number is obtained by averaging all estimated rates v and their propagated standard deviations σvv, but note that the size, distribution and representativeness of the observed glaciated areas are not taken into account. For this particular result, the absolute value of the estimated error, 0.47 m a−1 is larger than the estimate rate at −0.17 m a−1 . This result indicates that, given the measurements, it is most likely that glaciers on the Tibetan plateau were thinning between 2003 and 2009, but there is some significant chance that they were actually thickening. A more extensive study on the uncertainties associated to glacier mass balance studies from a geostatistics perspective can be found in [34].
The largest decrease in glacial thickness occurred at the Hengduan Mountains, compare Figure 5. The estimated rate equals −2.03 ± 0.73 m a−1 with a RMSE of 0.32 m. The observed glaciated area consists of two GLIMS glaciers facing east. Although there are little discrepancies between the GLIMS glacier outlines and the Landsat 8 OLI/TIRS, captured on 13 August 2013, Figure 5 indicates that glaciers have retreated significantly between ~2002, the time corresponding to the GLIMS database, and 2013. On the other hand, the observed glaciated area facing north at Western Kunlun Mountains had an elevation increase rate of 1.25 ± 0.51 m a−1 and a RMSE of 3.09 m, as illustrated in Figure 6. Overlaying the GLIMS glacier mask on the Landsat 8 OLI/TIRS image from 18 September 2013 indicates that in this area the glacier extent is relatively stable.
For each basin belonging to the Tibetan Plateau, a mean thinning or thickening rate v B ¯ ± σ ¯ B is estimated, as average of rates v and propagated standard deviations σvv. The result is shown in Table 3. In practice, the rate per basin is of course affected by the area of each glacier within the basin. However, in this study we only estimate trends representative of nearby glacier groups. A next but far from trivial step would be to design an interpolation scheme taking the sparsely available trends as input and use them to estimate an overall trend while incorporating e.g., the relative location, orientation, and representativeness of each available trend. Here, the area of glaciers is not taken into account when estimating overall glacial rates. The results show that mass loss due to glacier-thinning seems to take place in most of the basins, excluding Tarim Basin. Subsequently, lost or gained water volumes from glaciers by basin are approximately estimated, by multiplying the mean glacial thickness change rate with the total glacier area of each basin, as shown in Table 3.

3.2. Impact of Orientation on Glacial Thickness Change

The results indicate that glacial thickness change indeed strongly depends on the relative position in a mountain range. Most glaciers at a north face increase in volume, although some decrease but in that case at a slower rate than its south-facing counterpart. In total, there are 15 pairs of observed glaciated areas, i.e., adjacent glaciated areas located on opposite sides of the main mountain ridge, all listed in Table 4. Such situation is illustrated in Figure 7, showing the Western Kunlun Mountains range. The temporal trends between 2003 and 2009 on the north-facing glaciated area A equaled 0.69 ± 0.30 m a−1 while on its south-facing counterpart, glaciated area B, the trend had opposite sign, equaling −1.02 ± 0.29 m a−1. Similarly, the glacial thickness change rates at E, facing north, and F, facing southeast, were 0.58 ± 0.28 m a−1 and −0.29 ± 0.44 m a−1, respectively. Furthermore, the glacial thickness on C, toward the northeast, was estimated to decrease at a rate of 0.09 ± 0.30 m a−1 while glaciers in area D, toward the southwest, thinned at a rate of −0.29 ± 0.20 m a−1. A possible explanation is that south-facing glaciers receive much more solar radiation than north-facing glaciers. Even glaciated area C, oriented toward the northeast, faces the sun more than areas A and E. Similarly, glaciated area D, oriented toward the southwest, is receiving less sunlight than glaciated areas B and F. Additionally, this can be also the effect of precipitation driven by orography.

4. Discussion

In this section, we discuss the sensitivity of our results to the removal of ICESat footprints based on terrain surface criteria and the GLIMS glacier mask. First we discuss the impact of the terrain surface criteria for assessing the signal quality of the ICESat measurements. Second, the GLIMS glacier mask is static which has some effect on the estimation of glacial thickness change trend. Finally a comparison of our result to previous research is presented.

4.1. Exploring Terrain Surface Criteria

Several large glaciers sampled by ICESat footprints were considered to explore appropriate terrain surface criteria. The following relations were studied while determining the thresholds for terrain slope and roughness: glacier elevation difference Δh vs. slope S, roughness R and elevation hSRTM, respectively; and slope S vs. elevation hSRTM. The results are illustrated here for one case study considering a glacier area at the Gurla Mandhata I Mountains The results indicate that glacier elevation differences Δh increase with terrain slope, as illustrated in Figure 8a. The existence of such a slope bias is already described [35]. Large valley glaciers often have a surface roughness of below 20 m, see Figure 8b. In addition, a larger surface roughness will result in a positive bias in the estimated glacial thickness.
The relaxation of the slope threshold results in an increase in the number of accepted ICESat track segments sampling a glaciated area. This is illustrated in Figure 9 for an area in the Hengduan Mountains (No. 6 in Table S1). In Figure 9a, a number of 10 track segments was accepted, given a slope threshold of 15°. Based on these track segments, a trend was estimated with a RMSE of 4.18 m. In Figure 9b, the slope threshold was relaxed to 25°, resulting in a total number of 13 track segments. However, the quality of the final trend (RMSE = 6.39 m) decreases with the increase of the number of track segments. These two examples show some of the impacts of the slope and roughness thresholds.
In previous research, the results were annual glacial thickness change trends for defined regions [13,16]. These trends were directly estimated from all glacier elevation differences between ICESat elevations and the reference SRTM DEM on glacier areas, after removing footprints affected by clouds. This method ensures the availability of sufficient ICESat footprints to estimate trends in glacial thickness for relatively large regions. However, it ignores the impact of the high relief terrain characteristics of the Tibetan Plateau and surrounding mountain ranges. In addition, their definition of the sampled regions somehow smooths out significant signal, as it lumps together glaciers with different characteristics with respect to orography and orientation. Clearly there is a difficult trade-off between using more elevations of less individual quality against using less elevations of better quality.

4.2. State of the GLIMS Glacier Mask

Observations serving as input for the GLIMS glacier mask were obtained from 1978 to 2002, using aerial photographs, topographic maps and in situ measurements [24]. Because of remoteness and harsh climatic conditions on the Tibetan Plateau, it is difficult to make field investigation, therefore the Chinese glacier inventory that was used to establish the GLIMS glacier mask took place at different periods. The inventory was organized per drainage basin. The inventory for example took place at Qilian Mountains in 1981, at the Inner Plateau in 1988, etc. Positional uncertainty is expressed as a distance of 20 m, i.e., a given location lies within a circle of 20 m radius from the true location. In addition, recent studies report that the total glacier area on the Tibetan Plateau is shrinking [2,4,5,7,8,9]. Therefore, in this study some ICESat footprints acquired between 2003 and 2009 will fall within the GLIMS glacier outlines but are not sampling a real glacier anymore. This will affect the mean elevation difference Δ h ¯ i at the ICESat acquisition time ti. However, the number of such footprints within the same ICESat track segment is not large because the along track distance between consecutive footprints is approximately 170 m, and criteria on terrain surface are in place to remove uncertain footprints.
To further improve the glacial thickness change trends derived from ICESat/GLAS data, two techniques were applied. First, the glacier mask could be checked for each ICESat campaign using contemporary spectral (e.g., Landsat 8) or SAR data (e.g., Sentinel 1). Alternatively, classification techniques could be applied to the ICESat full waveform signals (GLA01 or GLA06 product) to verify if a ICESat signal is sampling snow, ice or rock [36]. Applying both types of analysis for the complete Tibetan Plateau is quite labor intensive however. Additionally, the most cloud free Landsat scenes, acquired between 2003 and 2011 to delineate glacier outlines [13,16]. However, it is difficult to match the acquisition time of ICESat campaigns with Landsat data for the full observed period for the whole Tibetan Plateau.

4.3. Glacial Thickness Changes for Sub-Regions

Our results consider annual glacial thickness change trends for relatively small areas. It is interesting to compare it with previous research. Neckel et al. (2014) grouped glaciers on the Tibetan Plateau into eight sub-regions, as illustrated in Figure 10 [16]. One of their results consists of annual glacial thickness change trends for each of these eight sub-regions. Accordingly we estimated glacial thickness change trends for the same eight sub-regions as well. For each sub-region, a mean glacial thickness change rate v ¯ R ± σ ¯ R is estimated as average of the glacial thickness change rates v and propagated standard deviations σvv of the observed glaciated areas within the sub-region. The results are presented in Table 5 and compared to Neckel’s Δh trends.
The comparison indicates that sub-regions (A, F, G, and H), relatively densely covered by glaciers, have a similar thickness change rate. Considering the other sub-regions, sub-region D has a somehow similar trend while rates in sub-regions B and C have a relative large disparity. The disparity between sub-regions B and C may be caused by: (i) the low number of observed glaciated areas; and (ii) differences in orientation of the observed glaciated areas: sub-region B consists of two south-facing glaciated areas and one north-facing glaciated area while sub-region C consists of three south-facing glaciated areas and two north-facing glaciated areas. At sub-region E, in case we set S0 = 20° and R0 = 15 m, the number of ICESat footprints is not enough to estimate a temporal trend. We assume that the total number of observed glaciated areas per sub-region and their orientation affect these mean glacial thickness change rates. That is, when the number of observed glaciated areas is large enough and observed glaciated areas located on opposite sides of the main mountain ridge are approximate, the mean glacial thickness change trend per sub-region is going to be more reliable.
Generally, our results are comparable to elevation change rates v ¯ G ± σ ¯ G estimated for high-mountain Asian glaciers by Gardner et al. (2013) [15]. Both results indicate that most of the glaciers in the Tibetan Plateau are thinning, except for the Western Kunlun Mountains, as shown in Table 6. The strongest glacier-thinning occurs in the Himalaya range and in the Hengduan mountains. The glacial thickness change rate in the western and inner plateau is near balanced or nearly equals zero. Inversely glaciers in the Western Kunlun Mountains are thickening.

4.4. Representativeness of An Observed Glaciated Area

A difficult question is to what extent the sparse estimates obtained by ICESat are representative for the full population of the Tibetan Plateau glaciers. This question cannot be answered here but we can assess which fraction of the glaciers is sampled. For this purpose, we determine the ratio κ between glaciated area sampled by ICESat footprints and the total glaciated area, following Equation (13). Here N is the total number of accepted ICESat footprints, AF is the area covered by one ICESat footprint and AG is the total sampled glaciated area.
κ = N × A F A G
A glaciated area can be considered to be well sampled if the total number of ICESat footprints sampling is large, while its total area is relatively small. An ICESat footprint with its diameter of 70 m occupies an area AF of ~3850 m2. For example, in Figure 2, glaciated area A occupies 30.6 km2 and is sampled by 108 accepted ICESat footprints. Therefore, A’s sample ratio equals 0.0136. Similarly, glaciated area B occupies 8.5 km2 and is sampled by 94 accepted ICESat footprints, so B’s sample ratio is 0.0426. In this way, the sample ratio for each of 90 observed glaciated areas is determined (see Table S1). Note that this ratio does not take the spatial and temporal distribution of the ICESat footprints into account, and therefore only provides a very rough indication on how well a glaciated area is sampled.
Similarly, the sample ratio for all observed glaciated areas on the whole Tibetan Plateau could be computed as well. As a result, the total area of 90 observed glaciated areas for the whole Tibetan Plateau is 5831.5 km2 and these glaciated areas were sampled by a total number of 16,002 accepted ICESat footprints. Thus in this case the sample ratio equals 0.0106. Note that one location might be sampled by several ICESat footprints from different epochs. That is not taken into account in this first assessment.

5. Conclusions

By exploiting ICESat laser altimetry data, thickness change rates of 90 glaciated areas on the whole Tibetan Plateau were estimated between 2003 and 2009. By considering terrain surface criteria slope and roughness, temporal glacial thickness change trends for the whole Tibetan Plateau were evaluated for 15 different settings. The results show that the settings of terrain slope and roughness equaling 20° and 15 m to remove uncertain ICESat footprints, respectively, are appropriate for the steep and rough glaciers in the Tibetan Plateau. In addition, the orientation of glaciers has been taken into account. The study indicated that most of the observed glaciated areas in the Himalaya, the Hengduan Mountains and the Tanggula Mountains experienced a serious thinning while in most of the observed areas in the Western Kunlun Mountains north-facing glaciers were thickening while south-facing glaciers were thinning. In addition, glacial thickness changes indeed strongly depend on the relative position in a mountain range. Most north-facing glaciers increase in thickness, although some decrease but, in those cases, at a slower rate than their south-facing counterpart.
Our results complement previously estimated water level changes of Tibetan lakes [37,38]. Using additional explicit runoff relations between glaciers and lakes [39], correlations between glacial and lake level changes can be determined to improve understanding of the water balance on the Tibetan Plateau.

Supplementary Materials

The following are available online at www.mdpi.com/2072-4292/9/2/160/s1, Table S1: Rates of glacial thickness changes on the Tibetan Plateau between 2003 and 2009.

Acknowledgments

This work was jointly supported by the EU-FP7 project CEOP-AEGIS (Grant No. 212921) and by the Vietnam Ministry of Education and Training.

Author Contributions

This work is a part of the PhD research performed by Vu Hien Phan, under the supervision of Massimo Menenti and Roderik Lindenbergh.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. National Snow and Ice Data Center. GLIMS Glacier Database, Version 1; NSIDC (National Snow and Ice Data Center): Boulder, CO, USA, 2005. Available online: http://0-dx-doi-org.brum.beds.ac.uk/10.7265/N5V98602 (accessed on 10 September 2013).
  2. Yao, T.; Thompson, L.; Yang, W.; Yu, W.; Gao, Y.; Gou, X.; Yang, X.; Duan, K.; Zhao, H.; Xu, B.; et al. Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings. Nat. Clim. Chang. 2012, 2, 663–667. [Google Scholar] [CrossRef]
  3. Sorg, A.; Bolch, T.; Stoffel, M.; Solomina, O.; Beniston, M. Climate change impacts on glaciers and runoff in Tien Shan (Central Asia). Nat. Clim. Chang. 2012, 2, 725–731. [Google Scholar] [CrossRef]
  4. Wang, P.; Li, Z.; Gao, W. Rapid shrinking of glaciers in the Middle Qilian Mountain region of Northwest China during the last 50 years. J. Earth Sci. 2011, 22, 539–548. [Google Scholar] [CrossRef]
  5. Tian, H.; Yang, T.; Liu, Q. Climate change and glacier area shrinkage in the Qilian mountains, China, from 1956 to 2010. Ann. Glaciol. 2014, 55, 187–197. [Google Scholar] [CrossRef]
  6. Bolch, T.; Yao, T.; Kang, S.; Buchroithner, M.F.; Scherer, D.; Maussion, F.; Huintjes, E.; Schneider, C. A glacier inventory for the western Nyainqentanglha Range and the Nam Co 15 Basin, Tibet, and glacier changes 1976–2009. Cryosphere 2010, 4, 419–433. [Google Scholar] [CrossRef]
  7. Wei, J.; Liu, S.; Gou, W.; Yao, X.; Xu, J.; Bao, W.; Jiang, Z. Surface-area changes of glaciers in the Tibetan Plateau interior area since the 1970s using recent Landsat images and historical maps. Ann. Glaciol. 2014, 55, 213–222. [Google Scholar] [CrossRef]
  8. Zhang, Y.; Liu, S.; Xu, J.; Shangguan, D. Glacier change and glacier runoff variation in the Tuotuo River Basin, the source region of Yangtze River in western China. Environ. Geol. 2008, 56, 59–68. [Google Scholar] [CrossRef]
  9. Ye, Q.; Zhong, Z.; Kang, S.; Stein, A.; Wei, Q.; Liu, J. Monitoring glacier and supra-glacier lakes from space in Mt. Qomolangma Region of the Himalayas on the Tibetan Plateau in China. J. Mt. Sci. 2009, 6, 211–220. [Google Scholar] [CrossRef]
  10. Gardelle, J.; Berthier, E.; Arnaud, Y. Slight mass change of Karakoram glaciers in the early twenty-first century. Nat. Geosci. 2012, 5, 322–325. [Google Scholar] [CrossRef]
  11. Wang, D.; Kääb, A. Modeling glacier elevation change from DTM time series. Remote Sens. 2015, 7, 10117–10142. [Google Scholar] [CrossRef]
  12. Levinsen, J.F.; Khvorostovsky, K.; Ticconi, F.; Shepherd, A.; Forsberg, R.; Sørensen, L.S.; Muir, A.; Pie, N.; Felikson, D.; Flament, T.; et al. ESA ice sheet CCI: Derivation of the optimal method for surface elevation change detection of the Greenland ice sheet–round robin results. Int. J. Remote Sens. 2015, 36, 551–573. [Google Scholar] [CrossRef]
  13. Kääb, A.; Berthier, E.; Nuth, C.; Gardelle, J.; Arnaud, Y. Contrasting patterns of early twenty-first-century glacier mass change in the Himalayas. Nature 2012, 488, 495–498. [Google Scholar] [CrossRef] [PubMed]
  14. Kropácek, J.; Neckel, N.; Bauder, A. Estimation of volume changes of mountain glaciers from ICESat data: An example from the Aletsch Glacier, Swiss Alps. Cryosphere Discuss. 2013, 7, 3261–3291. [Google Scholar] [CrossRef]
  15. Gardner, A.S.; Moholdt, G.; Cogley, J.G.; Wouters, B.; Arendt, A.A.; Wahr, J.; Berthier, E.; Hock, R.; Pfeffer, W.T.; Kaser, G.; et al. A reconciled estimate of glacier contributions to sea level rise: 2003 to 2009. Science 2013, 340, 852–857. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Neckel, N.; Kropacek, J.; Bolch, T.; Hochschild, V. Glacier mass changes on the Tibetan Plateau 2003–2009 derived from ICESat laser altimetry measurements. Environ. Res. Lett. 2014, 9, 014009. [Google Scholar] [CrossRef]
  17. Quincey, D.J.; Luckman, A.; Benn, D. Quantification of Everest region glacier velocities between 1992 and 2002, using satellite radar interferometry and feature tracking. J. Glaciol. 2009, 55, 596–606. [Google Scholar] [CrossRef]
  18. Zwally, H.J.; Schutz, R.; Bentley, C.; Bufton, J.; Herring, T.; Minster, J.; Spinhirne, J.; Thomas, R. GLAS/ICESat L2 Global Land Surface Altimetry Data, Version 31; NASA National Snow and Ice Data Center Distributed Active Archive Center: Boulder, CO, USA, 2010; Available online: http://0-dx-doi-org.brum.beds.ac.uk/10.5067/ICESAT/GLAS/DATA227 (accessed on 10 June 2010).
  19. Schutz, B.E. Laser Footprint Location (Geolocation) and Surface Profiles; Center for Space Research, The University of Texas: Austin, TX, USA, 2002. [Google Scholar]
  20. Duong, H.; Pfeifer, N.; Lindenbergh, R.C. Full waveform analysis: ICESat laser data for land cover classification. In Proceedings of the ISPRS Mid-term Symposium, Remote Sensing: From Pixel to Process, Enschede, The Netherlands, 8–11 May 2006; pp. 31–35.
  21. Jarvis, A.; Reuter, H.I.; Nelson, A.; Guevara, E. Hole-Filled SRTM for the Globe Version 4. 2008. The CGIAR-CSI SRTM 90m Database. Available online: http://srtm.csi.cgiar.org (accessed on 10 September 2013).
  22. Zandbergen, P. Applications of shuttle radar topography mission elevation data. Geogr. Compass 2008, 2, 1404–1431. [Google Scholar] [CrossRef]
  23. Mukul, M.; Srivastava, V.; Mukul, M. Accuracy analysis of the 2014–2015 Global Shuttle Radar Topography Mission (SRTM) 1 arc-sec C-Band height model using International Global Navigation Satellite System Service (IGS) Network. J. Earth Syst. Sci. 2015, 125, 909–917. [Google Scholar] [CrossRef]
  24. Frey, H.; Paul, F. On the suitability of the SRTM DEM and ASTER GDEM for the compilation of topographic parameters in glacier inventories. Int. J. Appl. Earth Obs. Geoinf. 2012, 18, 480–490. [Google Scholar] [CrossRef]
  25. Nuth, C.; Kääb, A. Co-registration and bias correction of satellite elevation data sets for quantifying glacier thickness change. Cryosphere 2011, 5, 271–290. [Google Scholar] [CrossRef]
  26. Shi, Y.; Liu, C.; Kang, E. The glacier inventory of China. Ann. Glaciol. 2009, 50, 1–4. [Google Scholar]
  27. Guo, W.; Liu, S.; Xu, J.; Wu, L.; Shangguan, D.; Yao, X.; Wei, J.; Bao, W.; Yu, P.; Liu, Q.; et al. The second Chinese glacier inventory: Data, methods and results. J. Glaciol 2015, 61, 357–372. [Google Scholar] [CrossRef]
  28. Verdin, K.L.; Godt, J.W.; Funk, C.; Pedreros, D.; Worstell, B.; Verdin, J. Development of A Global Slope Dataset for Estimation of Landslide Occurrence Resulting from Earthquakes; US Geological Survey: Colorado, CO, USA, 2007.
  29. Shi, J.; Menenti, M.; Lindenbergh, R.C. Parameterization of surface roughness based on ICESat/GLAS full waveforms: A case study on the Tibetan Plateau. J. Hydrometeorol. 2013, 14, 1278–1292. [Google Scholar] [CrossRef]
  30. Sinnott, R.W. Virtues of the haversine. Sky Telesc. 1984, 68, 159. [Google Scholar]
  31. Lay, D.C. Orthogonality and Least Squares. In Linear Algebra and Its Applications, 3rd ed.; Addison Wesley: Boston, FL, USA, 2002; pp. 329–392. [Google Scholar]
  32. Pavlis, N.K.; Holmes, S.A.; Kenyon, S.C.; Factor, J.K. An Earth Gravitational Model to 10 Degree 2160: EGM2008. In Proceedings of the 2008 General Assembly of the European Geosciences Union, Vienna, Austria, 13–18 April 2008.
  33. Teunissen, P.J.G. The model with observation equations. In Adjustment Theory: An Introduction; VSSD: Delft, The Netherlands, 2003; pp. 39–60. [Google Scholar]
  34. Rolstad, C.; Haug, T.; Denby, D. Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: Application to the western Svartisen ice cap. Nor. J. Geogr. 2009, 55, 666–680. [Google Scholar] [CrossRef]
  35. Slobbe, D.C.; Lindenbergh, R.C.; Ditmar, P. Estimation of volume change rates of Greenland’s ice sheet from ICESat data using overlapping footprints. Remote Sens. Environ. 2008, 112, 4204–4213. [Google Scholar] [CrossRef]
  36. Molijn, R.; Lindenbergh, R.C.; Gunter, B. ICESat laser full waveform analysis for the classification of land cover types over the cryosphere. Int. J. Remote Sens. 2011, 32, 8799–8822. [Google Scholar] [CrossRef]
  37. Zhang, G.Q.; Xie, H.J.; Kang, S.C.; Yi, D.H.; Ackley, S.F. Monitoring lake level changes on the Tibetan Plateau using ICESat altimetry data (2003–2009). Remote Sens. Environ. 2011, 115, 1733–1742. [Google Scholar] [CrossRef]
  38. Phan, V.H.; Lindenbergh, R.C.; Menenti, M. ICESat derived elevation changes of Tibetan lakes between 2003 and 2009. Int. J. Appl. Earth Obs. 2012, 17, 12–22. [Google Scholar] [CrossRef]
  39. Phan, V.H.; Lindenbergh, R.C.; Menenti, M. Geometric dependency of Tibetan lakes on glacial runoff. Hydrol. Earth Syst. Sci. 2013, 17, 4061–4077. [Google Scholar] [CrossRef]
Figure 1. GLIMS glacier outlines and ICESat L2D-campaign tracks superimposed on the SRTM DEM over the Tibetan Plateau.
Figure 1. GLIMS glacier outlines and ICESat L2D-campaign tracks superimposed on the SRTM DEM over the Tibetan Plateau.
Remotesensing 09 00160 g001
Figure 2. The ICESat-sampled glaciers having similar orientation were grouped into glaciated areas A and B in the Western Kunlun Mountains.
Figure 2. The ICESat-sampled glaciers having similar orientation were grouped into glaciated areas A and B in the Western Kunlun Mountains.
Remotesensing 09 00160 g002
Figure 3. Distributions of the mean elevation differences and temporal glacial thickness change trends between 2003 and 2009 at the glaciated areas A and B.
Figure 3. Distributions of the mean elevation differences and temporal glacial thickness change trends between 2003 and 2009 at the glaciated areas A and B.
Remotesensing 09 00160 g003
Figure 4. Glacial thickness change rates on the Tibetan Plateau between 2003 and 2009.
Figure 4. Glacial thickness change rates on the Tibetan Plateau between 2003 and 2009.
Remotesensing 09 00160 g004
Figure 5. The maximal rate of glacial thickness decrease between 2003 and 2009 at Mount Hengduan. The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 13 August 2013.
Figure 5. The maximal rate of glacial thickness decrease between 2003 and 2009 at Mount Hengduan. The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 13 August 2013.
Remotesensing 09 00160 g005
Figure 6. Strong glacial thickening between 2003 and 2009 at Western Kunlun Mountains The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 18 September 2013.
Figure 6. Strong glacial thickening between 2003 and 2009 at Western Kunlun Mountains The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 18 September 2013.
Remotesensing 09 00160 g006
Figure 7. Different rates of glacial thickness changes between 2003 and 2009 at the north and south face of the Western Kunlun Mountains The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 11 September 2013, and adding the locations of observed glaciated areas with thickness change rates.
Figure 7. Different rates of glacial thickness changes between 2003 and 2009 at the north and south face of the Western Kunlun Mountains The figure is created by overlaying the GLIMS glacier outlines on the Landsat 8 OLI/TIRS image from 11 September 2013, and adding the locations of observed glaciated areas with thickness change rates.
Remotesensing 09 00160 g007
Figure 8. Relations between: (a) glacier elevation difference and slope; and (b) glacier elevation difference and roughness. Glacier elevation differences are between ICESat campaigns L2A, L3A, L3D and L3G and the SRTM DEM reference surface over a glaciated area (No. 20 in Table S1) at Mount Guala Mandhata I, belonging to the Ganges Basin.
Figure 8. Relations between: (a) glacier elevation difference and slope; and (b) glacier elevation difference and roughness. Glacier elevation differences are between ICESat campaigns L2A, L3A, L3D and L3G and the SRTM DEM reference surface over a glaciated area (No. 20 in Table S1) at Mount Guala Mandhata I, belonging to the Ganges Basin.
Remotesensing 09 00160 g008
Figure 9. Estimations of glacial thickness change trends with varying slope S0 thresholds: (a) 15°; and (b) 25° at a glaciated area (No. 6 in Table S1) in the Hengduan Mountains, belonging to the Brahmaputra Basin. In this example the roughness R0 was kept fixed at 15 m.
Figure 9. Estimations of glacial thickness change trends with varying slope S0 thresholds: (a) 15°; and (b) 25° at a glaciated area (No. 6 in Table S1) in the Hengduan Mountains, belonging to the Brahmaputra Basin. In this example the roughness R0 was kept fixed at 15 m.
Remotesensing 09 00160 g009
Figure 10. Sub-regions applied for monitoring glacial thickness change at the Tibetan Plateau.
Figure 10. Sub-regions applied for monitoring glacial thickness change at the Tibetan Plateau.
Remotesensing 09 00160 g010
Table 1. The attributes related to each ICESat measurement.
Table 1. The attributes related to each ICESat measurement.
NameAttribute Description
TimeICESat acquisition time or arrival time of the laser pulse on the reflecting surface in UTC “dd-mm-yyyy” format, derived from the GLA14 data
LatGeodetic latitude in degrees, derived from the GLA14 data
LonGeodetic longitude in degrees, derived from the GLA14 data
ElevElevation in meters above WGS84, derived from the GLA14 data
GdHtGeoid height in meters in the EGM2008 datum, derived from the GLA14 data
SatFlgSaturation correction flag, identifying possible saturation issues, derived from the GLA14 data
NumPkNumber of peaks in the Gauss waveform decomposition from the return echo, derived from the GLA14 data
SRTM_elevElevation in meters above EGM1996, derived from the SRTM DEM
SSurface slope in degrees, derived from the SRTM DEM
RSurface roughness in meters, derived from the SRTM DEM
96_08_HtGeoid height difference between EGM96 and EGM2008 in meters, derived from the NGA/NASA Earth Gravitational Model tools [ 32]
GIDIdentification code of the observed glaciated area
Table 2. Settings of terrain surface parameters for filtering ICESat footprints.
Table 2. Settings of terrain surface parameters for filtering ICESat footprints.
SettingS0 (Degree)R0 (m)N v ¯ (m a−1) σ ¯ v v (m a−1) R M S E ¯ (m)N5
S110533−0.210.202.9329
S215538−0.230.213.2634
S320543−0.120.213.0640
S4255490.010.233.3443
S5305540.040.234.0041
S6101037−0.250.252.8533
S7151055−0.060.332.9949
S8201076−0.020.393.7062
S92510980.130.444.2968
S103010117−0.040.455.4067
S11101539−0.210.262.8936
S12151563−0.150.403.0558
S132015900.170.474.0267
S142515122−0.210.564.8964
S153015146−0.210.615.9257
Table 3. Mean glacial thickness change rate per basin, where N is the number of observed glaciated areas and the total glacier area is obtained from the GLIMS glacier mask. Lost or gained water volumes from glaciers are approximately estimated, by multiplying the mean glacial thickness change rate with the total glacier area of each basin.
Table 3. Mean glacial thickness change rate per basin, where N is the number of observed glaciated areas and the total glacier area is obtained from the GLIMS glacier mask. Lost or gained water volumes from glaciers are approximately estimated, by multiplying the mean glacial thickness change rate with the total glacier area of each basin.
BasinTotal Glacier Area (km2)N v ¯ B ± σ ¯ B (m a−1)Water Volume (Gt a−1)
Brahmaputra16,0199−0.56 ± 0.49−8.97 ± 7.79
Ganges40338−0.99 ± 0.47−4.01 ± 1.90
Indus24095−0.03 ± 0.34−0.08 ± 0.82
Inner plateau870223−0.16 ± 0.48−1.39 ± 4.14
Salween18511−0.78 ± 0.81−1.44 ± 1.51
Tarim20,996390.21 ± 0.474.31 ± 9.79
Yangtze20125−1.14 ± 0.46−2.30 ± 0.93
Total56,56190−0.17 ± 0.47−9.62 ± 26.41
Table 4. List of pairs of glaciated areas that are adjacent, but located on opposite sides of the main mountain ridge. Here, Nf is the total number of accepted footprints. Locations A, B, C, D, E and F are indicated in Figure 7.
Table 4. List of pairs of glaciated areas that are adjacent, but located on opposite sides of the main mountain ridge. Here, Nf is the total number of accepted footprints. Locations A, B, C, D, E and F are indicated in Figure 7.
LatitudeLongitudeBasinOri.Nf v ± σ v v (m a−1)RMSE (m)
128.18490.544BrahmaputraS261−0.09 ± 0.398.68
228.24890.543BrahmaputraN71−0.14 ± 0.407.13
328.26186.296GangesS323−1.83 ± 0.373.40
428.33686.302GangesN930.12 ± 0.254.64
530.41581.306GangesS80−0.90 ± 0.695.83
630.46981.310GangesN99−0.74 ± 0.543.40
730.93683.494Inner plateauE831.63 ± 0.589.21
831.02283.468Inner plateauW160−0.46 ± 0.363.56
933.91390.659Inner plateauS92−0.47 ± 0.203.92
1033.95490.670YangtzeN342−0.60 ± 0.303.23
1134.02479.763IndusSW79−1.38 ± 0.432.73
1234.05379.788IndusE185−0.07 ± 0.201.51
1334.28881.946Inner plateauS1061.23 ± 0.502.76
1434.32781.946Inner plateauN1680.21 ± 0.472.25
1535.28480.685Inner plateau (B)S998−1.02 ± 0.294.19
1635.52380.713Tarim (A)N13200.69 ± 0.303.38
1735.30181.430Inner plateau (D)SW635−0.29 ± 0.201.73
1835.38881.397Tarim (C)NE633−0.09 ± 0.301.44
1935.41081.612Tarim (F)SE338−0.44 ± 0.443.46
2035.50881.624Tarim (E)N3800.58 ± 0.281.79
2135.47082.143Inner plateauS92−1.50 ± 0.794.41
2235.51682.162Tarim N77−1.02 ± 0.435.07
2335.65585.620Inner plateauS1181.82 ± 0.485.08
2435.69685.613Inner plateauN257−0.04 ± 0.242.85
2535.77477.130Tarim W930.06 ± 0.574.74
2635.81277.148Tarim N470.19 ± 0.573.16
2736.02490.962Tarim S428−0.80 ± 0.387.03
2836.09990.936Inner plateauN494−0.55 ± 0.222.88
2936.77384.903Inner plateauS59−0.13 ± 0.562.89
3036.81384.895Tarim N520.03 ± 0.782.44
Table 5. Mean glacial thickness change rates per sub-region, where N is the number of observed glaciated areas within each sub-region, compared to the results of Neckle et al. (2014) [16].
Table 5. Mean glacial thickness change rates per sub-region, where N is the number of observed glaciated areas within each sub-region, compared to the results of Neckle et al. (2014) [16].
Sub-RegionNameN v ¯ R ± σ ¯ R (m a−1)Neckel’s v ¯ G ± σ ¯ G (m a−1)
AWestern Kunlun Mountains200.16 ± 0.440.04 ± 0.29
BZangser Kangri and Songzhi Peak30.86 ± 0.310.44 ± 0.26
CQilian Mountains and Eastern Kunlun Mountains50.03 ± 0.47−0.90 ± 0.28
DTanggula Mountains and Dongkemadi Ice Cap6−0.88 ± 0.41−0.68 ±0.29
EWestern Nyainqentanglha range0NA−0.23 ± 0.33
FGangdise Mountains8−0.60 ± 0.50−0.44 ± 0.26
GCentral and Eastern Tibetan Himalaya8−0.70 ± 0.46−0.78 ± 0.27
HEastern Nyainqentanglha and Hengduan Mountains6−0.67 ± 0.58−0.81 ± 0.32
Table 6. Mean glacial thickness change rates per mountain region on the Tibetan Plateau, compared to the results of Gardner et al. (2013) [15].
Table 6. Mean glacial thickness change rates per mountain region on the Tibetan Plateau, compared to the results of Gardner et al. (2013) [15].
High Mountain Regions v ¯ R ± σ ¯ R (m a−1) Gardner’s v ¯ G ± σ ¯ G (m a−1)
The Himalaya range−0.81 ± 0.46
Western −0.53 ± 0.13
Central −0.44 ± 0.20
Eastern −0.89 ± 0.13
The Hengduan mountains−0.67 ± 0.58−0.40 ± 0.41
The western and inner plateau−0.05 ± 0.450.02 ± 0.14
The Western Kunlun Mountains0.20 ± 0.450.17 ± 0.15

Share and Cite

MDPI and ACS Style

Phan, V.H.; Lindenbergh, R.; Menenti, M. Assessing Orographic Variability in Glacial Thickness Changes at the Tibetan Plateau Using ICESat Laser Altimetry. Remote Sens. 2017, 9, 160. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9020160

AMA Style

Phan VH, Lindenbergh R, Menenti M. Assessing Orographic Variability in Glacial Thickness Changes at the Tibetan Plateau Using ICESat Laser Altimetry. Remote Sensing. 2017; 9(2):160. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9020160

Chicago/Turabian Style

Phan, Vu Hien, Roderik Lindenbergh, and Massimo Menenti. 2017. "Assessing Orographic Variability in Glacial Thickness Changes at the Tibetan Plateau Using ICESat Laser Altimetry" Remote Sensing 9, no. 2: 160. https://0-doi-org.brum.beds.ac.uk/10.3390/rs9020160

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

Article Metrics

Back to TopTop