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

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


Introduction
The Tibetan Plateau has steep and rough terrain and contains~37,000 glaciers, occupying an area of~56,560 km 2 [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 km 2 per individual glacier, corresponding to a mean rate of 2127 m 2 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.

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.

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.

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.

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.

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.

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.

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

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

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.

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. 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].
where the h i 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].
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ê i between the grid heights and the local 3 × 3 plane, best fitting in the least squares sense [31], following Equation (5).

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.

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

Time
ICESat acquisition time or arrival time of the laser pulse on the reflecting surface in UTC "dd-mm-yyyy" format, derived from the GLA14 data Lat Geodetic latitude in degrees, derived from the GLA14 data Lon Geodetic longitude in degrees, derived from the GLA14 data Elev Elevation in meters above WGS84, derived from the GLA14 data GdHt Geoid height in meters in the EGM2008 datum, derived from the GLA14 data SatFlg Saturation correction flag, identifying possible saturation issues, derived from the GLA14 data

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.
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. Here, we analyze different settings incorporating the terrain surface characteristics slope and roughness. We remove footprints with a slope S bigger than a threshold S 0 and roughness R bigger than a threshold R 0 . 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.

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 t i . The mean difference ∆h i and its standard deviation s i 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 t i and ∆h ij is the jth elevation difference, j = 1 ÷ k.
Each ICESat acquisition time t i 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 s i is less than a threshold Std 0 and the number of ICESat footprints Remote Sens. 2017, 9, 160 8 of 19 k is at least six footprints. The threshold Std 0 is defined to be equal to the roughness threshold R 0 for each setting with respect to terrain slope and roughness. To remove unreliable elevation differences, we build an iterative algorithm. That is, if s i is bigger than Std 0 and ∆h ij − ∆h i is maximal for j in 1 ÷ k, the jth elevation difference ∆h ij is removed. Then, ∆h i and s i are re-computed. This process is repeated until s i drops below Std 0 or k is less than six. In Figure 3, the values ∆h i and s i 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 S 0 , R 0 , and Std 0 are 15 • , 10 m, and 10 m, respectively.
where, = ∆ℎ ∆ℎ … ∆ℎ : the vector of the average elevation differences per epoch. = : the vector of parameters of the linear trend, offset x0 and rate v.
: 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 ̂= − . This value consists of a combination of possible data errors and mainly the non-validity of the linear regression model.
In addition, the propagated standard deviation σvv of the estimated rate v is given in Equation (11), where yy Q 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.

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]. where, : the vector of the average elevation differences per epoch.
x = x 0 v : the vector of parameters of the linear trend, offset x 0 and rate v.
: the design matrix, in which t i 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ê = y − Ax. This value consists of a combination of possible data errors and mainly the non-validity of the linear regression model.
Remote Sens. 2017, 9, 160 9 of 19 In addition, the propagated standard deviation σ vv of the estimated rate v is given in Equation (11), where Q yy denotes the variance matrix, in which s i is the standard deviation of the ith average difference. These values are considered as the confidence interval for the estimated glacial thickness change.
Continuing to the example of Figure

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 S 0 and roughness R 0 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 R 0 = 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, S 0 and R 0 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 σ vv are the resulting overall rate and its propagated standard deviation of glacial thickness change while RMSE is the average of the root mean square errors (RMSEs) of the linear regression model. Additionally, N 5 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 S 0 and R 0 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.

Overall Glacial Thickness Changes: Tibetan Plateau and Its Basins
In the case the thresholds S 0 = 20 • for terrain slope and R 0 = 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].

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

Basin
Total Glacier Area (km 2 ) N v B ± oe B (m a −1 ) Water Volume (Gt a −1 )

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.

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.

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 h SRTM , respectively; and slope S vs. elevation h SRTM . 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.

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.

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.  Table S1) at Mount Guala Mandhata I, belonging to the Ganges Basin.
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  Table S1) at Mount Guala Mandhata I, belonging to the Ganges Basin.
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.
Remote Sens. 2017, 9,160 15 of 20 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.  Table S1) in the Hengduan Mountains, belonging to the Brahmaputra Basin. In this example the roughness R0 was kept fixed at 15 m.

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 i h Δ 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.

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  Table S1) in the Hengduan Mountains, belonging to the Brahmaputra Basin. In this example the roughness R 0 was kept fixed at 15 m.

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 t i . 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.

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. 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.   [16].

Sub-Region
Name 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.   [16].

Sub-Region
Name 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 S 0 = 20 • and R 0 = 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.

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, A F is the area covered by one ICESat footprint and A G is the total sampled glaciated area.
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 m 2 . For example, in Figure 2, glaciated area A occupies 30.6 km 2 and is sampled by 108 accepted ICESat footprints. Therefore, A's sample ratio equals 0.0136. Similarly, glaciated area B occupies 8.5 km 2 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 km 2 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.

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.