High Mountain Asia glacier elevation trends 2003–2008, lake volume changes 1990–2015, and their relation to precipitation changes

We present an updated, spatially resolved estimate of 2003–2008 glacier volume changes for entire High Mountain Asia (HMA) from ICESat laser altimetry data. The results reveal a diverse pattern that is driven by spatially greatly varying glacier sensitivity, in particular to precipitation availability and changes. We introduce a spatially resolved zonation where ICESat samples are grouped into units of similar glacier behaviour, glacier type, and topographic settings. In several regions, our new zonation reveals local differences and anomalies that have not been described previously. A step-increase in precipi5 tation around 1997–2000 on the Tibetan Plateau (TP) caused thickening of glaciers in the Eastern Pamirs, Kunlun Shan and central TP by 0.1–0.7 m a−1. The thickening anomaly has a crisp boundary in the Eastern Pamir that continues just north of the central Karakoram. Glaciers in the south and east of the TP were thinning, with increasing rates towards southeast. The precipitation increase is reflected by growth of endorheic lakes in particular in the northern and eastern TP. We estimate lake volume changes through a combination of repeat lake extents from Landsat data and shoreline elevations from ICESat and 10 the SRTM DEM for over 1300 lakes. The rise in water volume contained in the lakes corresponds to 4–25 mm a−1, when distributed over entire catchments, for the areas where we see glacier thickening. The precipitation increase is also visible in sparse in-situ measurements and MERRA-2 climate reanalysis data, but less well in ERA Interim reanalysis data. Considering evaporation loss, the difference between average annual precipitation during the 1990s and 2000s suggested by these datasets is 34–100 mm a−1, depending on region, which can fully explain both lake growth, and glacier thickening (Kunlun Shan) or 15 glacier geometry changes (eastern TP). The precipitation increase reflected in these glacier changes possibly extended to the northern slopes of the Tarim Basin, where glaciers were nearly in balance in 2003–2008. Along the entire Himalaya, glaciers on the first orographic ridge, which are exposed to abundant precipitation, are thinning less than glaciers in the dryer climate of the inner ranges. Thinning rates in the Tien Shan vary spatially but are rather stronger than in other parts of HMA.

i m a l a y a n R a n g e Q i l i a n S h a n N y a i n q ê n ta nglha N y a i n q ê n ta nglha H e n g d u a n S .
T i e n S h a n B o r o h o r o T a n g g u l a T i b e t a n P l a t e a u P am ir A la i  Yao et al., 2012;Sakai et al., 2015). From the eastern Himalaya and southern/eastern TP to the northwest of HMA, there is a transition from predominant spring/summer accumulation to winter accumulation in the Hindu Kush and the western parts of Tien Shan (Palazzi et al., 2013;Bookhagen and Burbank, 2010;Rasmussen, 2013). Mountains in between, such as the Karakoram and western Himalaya, receive moisture from both sources (Kuhle, 1990;Bolch et al., 2012). The Kunlun Mountains, on the other hand, receive most precipitation around May . The seasonal timing of glacier 5 accumulation likely plays an important role for glacier sensitivity to a warming climate (Fujita, 2008;Mölg et al., 2012;Sakai and Fujita, 2017). Another important factor is total precipitation availability, which depends on continentality (Shi and Liu, 2000;Kuhle, 1990) and, on smaller spatial scale, on glacier location on or behind a mountain range that acts as a primary  Wagnon et al. (2013) and Sherpa et al. (2016) found indications of steep horizontal precipitation gradients within only a few kilometres on the outermost ridge of the Great Himalaya in the Khumbu region in Nepal. Vertical precipitation gradients at high altitude are still poorly understood. It is suggested that precipitation increases from dry mountain valley 5 bottoms to an elevation of 4000-6000 m a.s.l. and subsequently decreases again at even higher elevations (e.g. Immerzeel et al., 2014Immerzeel et al., , 2015. Many glaciers in HMA are debris-covered in their ablation areas, and the percentage of debris-covered ice varies greatly between different regions (Scherler et al., 2011;Gardelle et al., 2013). Recent studies have found that although debris-covered glaciers in HMA have stable front positions (Scherler et al., 2011), they melt on average just as fast as clean ice glaciers 10 Gardelle et al., 2012b;Pellicciotti et al., 2015). In this study, we distinguish thus not explicitly between debris-covered and debris-free glacier tongues.

Data and Methods
In this section we give a short overview of the data and methods used. Details can be found in the Appendix. As reference DEM for our ICESat processing and to derive lake shoreline elevations we use the DEM from the Shuttle Radar 20 Topography Mission. We used the C-band, non-void-filled SRTM DEM version at 3 arc-seconds resolution (SRTM3). As an alternative elevation reference, we used also the SRTM DEM at 1-arc-second resolution (SRTM1). We did not explore or use the recently published TanDEM-X global DEM as it was not yet available during our processing. Due to temporal inconsistency and substantial voids, we did also not use the ALOS PRISM World DEM (AW3D) or the WorldView satellite optical stereo HMA DEM. (Appendix A2).

25
As an estimate for regional and temporal precipitation patterns for the years 1980-2015 we use data from the Modern-Era Retrospective analysis for Research and Applications, version 2 (MERRA-2), available from the NASA Goddard Earth Sciences Data and Information Services Center, and ERA Interim at T255 spectral resolution. We use monthly summarised values of the variables total precipitation (PRECTOT / tp), snowfall (PRECSNO / sf) and evaporation (EVAP / e) from MERRA-2's surface flux diagnostics dataset tavg1_2d_flx_Nx and ERA Interim's Monthly Means of Daily Forecast Accumulations, 30 respectively. The High Asia Reanalysis (HAR), a product optimised for the TP region and with much finer spatial resolution, is unfortunately only available for the time period 2001-2011 which is too short for our study with respect to the lake volume changes investigated. Further, we use in-situ data from the five westernmost meteorological stations on the TP and Kunlun Shan 5 The Cryosphere Discuss., https://doi. org/10.5194/tc-2018-238 Manuscript under review for journal The Cryosphere Discussion started: 12 November 2018 c Author(s) 2018. CC BY 4.0 License.
( Fig. 1), provided by the China Meteorological Science Data Sharing Service Network. The data includes daily measurements of precipitation, mean air temperature, and for the four stations on the southwestern TP also evaporation. (Appendix A3).
We extract repeat lake coverage from the Global Surface Water dataset that is a classification of the entire Landsat archive into monthly and annual maps of surface water. The data is available within Google Earth Engine. Coverage is nearly complete 5 (>98%) starting from 2000 but considerably worse for some years of the 90s. (Appendix A4).

Methods for glacier volume change
We follow the double-differencing method explained in further detail in Kääb et al. (2012) and Treichler and Kääb (2016), with special consideration of issues mentioned in the above introduction. The difference between ICESat and SRTM elevations is further referred to as dh. Double differencing, i.e. fitting a linear trend through dh from several years, reveals how much method described in Treichler and Kääb (2016). Kääb (2016, 2017) found that ICESat clearly records the onset of winter snowfall in Norway during the split autumn 2008 campaign (stopped half way in mid-October and completed only in December). Analogue to Treichler and Kääb (2016), we estimate December 2008 snow bias from a linear regression of 5 October/December 2008 land dh on elevation and time. (Appendix B3)

Methods for lake volume change
In order to relate glacier changes and precipitation changes on the Tibetan Plateau to each other, and in particular to investigate if precipitation increases could be a reason for the positive glacier mass balances found in parts of the region, we also derive the volume changes of endorheic lakes on the TP. On first order, and by neglecting changes in subsurface water transport, additional 10 lake water masses should for endorheic lake systems either stem from increased lake inflow (mainly increased precipitation or enhanced glacier melt) or from reduced water loss (mainly changes in evaporation). Lake volume changes on the TP serve thus as potential proxies for precipitation changes, but help also to correct satellite gravimetric signals of glacier mass changes (see introduction; Appendix C).
We compute annual water volume change of the Tibetan lakes by multiplying annual lake areas with water level changes 15 from repeat water surface elevations for each year over the period 1990-2015. Maximum annual lake extents are obtained directly from the Global Surface Water data set. We retrieve the corresponding lake surface elevations in two ways: a) from SRTM DEM elevations of the lake shore by computing the median of interpolated DEM elevations for lake shore cells for each areal extent, and b) directly from ICESat footprint elevations on the lake areas for those lakes where ICESat data is available. To extend the lake elevation time series from method b) beyond the ICESat period of 2003-2009, we compute the 20 area-surface-elevation relationship for each lake by robust linear regression and apply this function to the areal extends of the years before and after the ICESat period. The so-extrapolated surface elevation values generate complete 1990-2015 time series for both areal extent and lake levels from SRTM and ICESat data, respectively. Our method is in parts similar to the methods used by previous studies but the inclusion of a DEM for deriving shoreline elevations, and thus lake water levels, in addition to altimetry data enabled us to produce volume change time series for one order of magnitude more lakes (>1300) 25 than derived previously. (Appendix C).
To minimise the effect of uncertainties in or erroneous estimates for individual years, we analyse time series in a summarised way through regression over time and as decadal averages, and apply a range of filters. (Appendix C1).
To estimate the lake water volume change in a way that can be related to glacier mass balances and precipitation changes (i.e. mm w.e. per m 2 ), we summarise and spatially distribute the water volume changes of all lakes within spatially confined 30 basins. These basins are based on endorheic catchments, but because many catchments only contain a single lake and exact catchment areas are not well defined on the TP (e.g., in very flat areas), we manually controlled, adjusted and aggregated the endorheic catchments of the USGS HydroSHEDS dataset at 15 arcsec spatial resolution to larger basins of comparable size and consisting of in average 5 catchments. (Appendix C2).

Methods for precipitation change
A change in precipitation, minus the part that is lost through evaporation and when neglecting changes in subsurface water 5 transport, should yield numbers that are directly scalable in relation to endorheic catchment water volume change or glacier mass balance, especially where the latter is governed by precipitation rather than temperature/melt. However, reanalysis data may not be very accurate in HMA due to a lack of ground measurements, and the few meteorological stations are not necessarily representative for a larger area. We therefore use raw precipitation data mainly to detect/confirm temporal and large-scale spatial patterns, and in a summarised way through decadal averages, rather than relying on annual numbers.  Figure 2a shows the 100 spatial units of glacier surface elevation change that result from the iterative manual zone delineation process. Spatial units needed to be large on the TP where glacier density is low, and could be rather small in the Karakoram which is intensely glacierised. Along major ridges such as the Himalaya, the units were designed narrow and along ridge 15 orientation in order to group glaciers under similar temperature and precipitation regimes rather than across orographic barriers.

Glacier thinning and thickening
Surface elevation change for the new spatial units and the 2°× 2°grid in Fig. 2b    20 out of 100 units we were not able to compute the potential biasing effect of December 2008 snow cover (e.g., due to lack of off-glacier samples). To ensure a consistent approach, we did therefore not apply this correction to the results presented above but instead added the difference due to bias correction to the error budget.

Lake changes on the TP
We receive valid (according to our filter procedures) water volume change time series for 89% of the median lake area (74% of all endorheic lakes) on the TP: 1009 lakes with SRTM-based lake surface elevations, thereof 103 also having ICESat-based 5 lake surface elevations (59% of the lake area). Extrapolated lake levels based on annual or campaign ICESat data yield the same results, but ICESat-based lake level change is on average 1.55 times larger than SRTM-based values. Multiplied with areal changes to receive volume changes, the relative difference is reduced to 1.09 times. Average 1990-2015 lake growth corresponds to 0.14 m a −1 (SRTM) and 0.18 m a −1 (ICESat) in lake-level change per year (Fig. 4a, robust linear regression of dV scaled with median lake area for easier comparison of values between lakes of very different size). All, except a handful 10 of lakes predominantly in the very south of the TP, grew during the studied time period, and growth of individual lakes is largest in the northern and eastern part of the TP. Figure 5 shows relative lake volume growth (based on SRTM lake levels) for individual lakes and regional medians over time for six regions: southwestern, eastern, central, northeastern, northwestern TP and Qilian Shan, indicated in Fig. 4a. (Note that the y-axis in Fig. 4a is relative to the total volume change dV over the time period observed and does not show absolute lake volumes; these are unknown as our method only detects lake level changes 15 and not lake depths). Rather than growing steadily, most lakes seem to have undergone a phase of sudden and rapid growth starting in ∼1997 and gradually slowing down until ∼2009, with rather stable conditions before and after this period. Relative lake volume change was most sudden and rapid for the northeastern, northwestern, central and eastern TP (the former three corresponding to areas with 2003-2008 glacier thickening). Lakes in the southern and southwestern part of the TP showed more varying and overall less growth, with a tendency to decrease after 2010. Endorheic lakes in the Qaidam basin/Qilian Shan 20 region further northeast also show a different and more varying evolution with slower growth that started only around 2004, but continued until ∼2012. The latter effect is also visible for the adjacent lakes on the northeastern TP (east Kunlun Shan). of the lake volume growth on the eastern TP due to considerably larger lake areas and lake density compared to the mostly small lakes further north/west. Table 1 shows additional water volumes accumulated between the two decades for the same  regions as above.
To yield values comparable to precipitation changes, the reader has to divide the total decadal differences dV given in the table by the number of years during which the additional water was accumulated. Assuming the change happened rather gradual during the entire decade, the specific annual water change would correspond to 1/10 of the values in Table   30 1. For instance, for water volumes using SRTM-based lake levels: 25±3 mm a −1 for the eastern TP, 4±1 mm a −1 for the southwestern TP, 6-7±1 mm a −1 for the central and northern TP, and 0.1±0.5 mm a −1 for the Qaidam basin/Qilian Shan region. Notably, there are considerable differences between catchments within each region (range for SRTM-based estimates: -5±1 to +35±6 mm a −1 , excluding one outlier of 163±7 mm a −1 for the catchment centred at 34.3°N / 88.8°E). The estimates based on SRTM-and ICESat lake levels aggregated for the six regions nevertheless agree very closely. The above annual values have to be doubled, or the dV values given in the table multiplied by 1/5, for instance, if one prefers to assume that the water volume increase happened during 5 years only, with stable conditions before and after -an assumption which seems also well plausible from Fig. 5.

5
Annual precipitation sums on the TP from meteorological stations range from as little as 50-100 mm a −1 (Shiqanhe and Thashkurgan stations, southwest TP and West Kunlun Shan) to 500-900 mm a −1 (Nielaer station, southern TP). Reanalysis values of both products used, MERRA-2 and ERA Interim, lie in between. All datasets record the majority of precipitation (>70%) during the monsoon-influenced summer months (May-September), except for Pulan and Nielaer, the two southernmost stations close to the Himalaya (only ca. 50% precipitation in summer).

10
Of the five meteorological stations available, especially Shiquanhe and Pulan show little change in precipitation and pan evaporation (Fig. 6). The Gaize station, located most central on the TP but still more south than our corresponding glacier unit, indicates a step-like precipitation increase around the year 2000, but data from only one station need of course to be interpreted with care due to potential local effects and changes to the station. A more gradual increase is visible in the Tashkurgan data.
Differences in decadal average precipitation range from −1 (Shiqanhe) to 60 mm (Nielaer) within 10 years, notably with 15 greatest relative change for the Gaize station (+42 mm per decade, a 25% increase) and Tashkurgan station (+16 mm or +22% per decade). Decadal differences are mostly (Nielaer, Tashkurgan) or exclusively (other stations) caused by an increase of precipitation during summer months. Pan evaporation reaches twice to tenfold of precipitation sums.
The two reanalyses used here differ considerably both in precipitation evolution and in estimated evaporation. Figures 8a (ERA Interim) and 8b (MERRA-2) show regional averaged annual sums for total precipitation, evaporation, and the difference of the two, for grid points within the TP lake catchment regions defined above. Notably, ERA Interim suggests considerably higher evaporation values than MERRA-2, in particular for the southwestern TP and the three northern regions, resulting in 5 much lower suggested net water availability in the areas where we see glacier thickening than it is the case for MERRA-2. Both reanalyses show an increase in precipitation starting from ca. 1995, but for ERA Interim, the evolution only lasts until ca. 2000 after which precipitation sums decrease. Also, the short-term precipitation increase is not visible for the northern parts of the TP, and it does not result in a noticeable decadal difference between average annual precipitation in 1990-1999 vs. 2000-2009. MERRA-2, on the other hand, rather suggests a step-wise increase with continuously higher precipitation sums until ca.
2010 for the entire TP, and even a continuous increase through 2015 for the northern part of the TP. For all six regions, this results in a total increase in precipitation of 34 mm (northwestern TP) to 100 mm (eastern TP) mm within 10 years. Except for 5 the Qilian Shan region, the change is exclusively driven by increasing summer precipitation (Fig. 7b). Winter precipitaion did not change noticeably (−9 to −2 mm decadal change for the five TP regions, +8 mm for Qilian Shan). When correcting these values with MERRA-2 estimates of actual evaporation, the total decadal increase in water availability ranges from 6 mm (central TP) to 68 mm (northeastern TP). The evaporation-corrected increase is even greater when looking at summer months only (31-77 mm per decade, compared to a decrease in water availability during winter months of −27 to 10 −6 mm, not shown). MERRA-2 suggests that 30-60% of precipitation on the TP falls as snow during the summer months and that the proportion of snow fall did not change noticeably between the decades (not shown).
The regions where MERRA-2 indicates increased summer precipitation correspond well with those areas on the TP and in Eastern Kunlun Shan with moderately negative to positive surface elevation change and/or endorheic lake growth. ERA Interim data indicates a similar pattern but with much lower (TP) precipitation increase, or even decrease in case of the Kunlun Shan region (Fig. 7a) Our above results on precipitation changes relate to decadal means in order to enable systematic comparison to other data. It is however important to note that these results vary if other time periods are chosen for aggregation. , for instance, summarize total annual precipitation amounts estimated from ERA-interim reanalysis over the Aru region, nortwestern TP, over 1979TP, over -1995TP, over and 1995TP, over -2008 to suggest a 33% increase between the latter both periods. where MERRA-2 data suggests a step-like increase of summer precipitation around the year 2000. The change in available precipitation amounts, lake water volume, and glacier mass balances are of the same magnitude and match well in terms of timing. The fact that glacier volumes are predominantly increasing in regions where also lake volumes increase, and the fact that lake volumes are also increasing in little or not glacierized basins, both suggest that the increases in lake volumes over the study region are not mainly driven by increased water influx from glacier mass loss (see e.g. Song et al., 2015). Though, glacier mass loss can certainly play an additional role for lakes with declining glaciers in their catchment. This is in line with, and extents geographically, water balance modelling by Lei et al. (2013) for six selected lakes in our East TP zone (Fig. 4b) 5 that suggests mainly precipitation increases to be behind the increases of lake volumes, accompanied by decreases in potential evaporation due to decreasing wind speed, and to a lesser extent increase in glacier runoff . For the Siling Co lake, potential evaporation, though decreasing overall over 1961-2010, showed actually stable conditions or a slight increase between the mid/end 1990s to 2010 (Guo et al., 2019), underlining even more the key role of precipitation increases for the observed lake volume increase.

10
On a local scale, and in contrast to the above regional view, there are considerable differences to previous findings in glacier changes, including the ones based on the same ICESat data. Compared to a visualisation of our results in a regular grid, we find that spatial aggregation matters: even within our study, only the manual zonation brings forward finer spatial differences e.g. from topographic-orographic setting. Our results also suggest that inconsistent sampling hypsometry, snow cover, and local vertical biases and elevation inconsistencies can have a severe biasing effect on ICESat-based glacier changes when not accounted for properly -in particular where they vary for different ICESat campaigns.

Precipitation increase on the TP
The spatial patterns of MERRA-2 decadal precipitation increases and the glacier growth on the TP and in the Kunlun Shan Both patterns correspond to the zonal boundary of 'extreme continental (polar) glaciers' suggested by Shi and Liu (2000), which encompasses the northwestern half of the TP, glaciers north of central Karakoram, the easternmost Pamir, and the entire 5 Kunlun Shan. On a coarser spatial and longer temporal scale, Kapnick et al. (2014) suggest that glacier accumulation in the Karakoram is least sensitive to atmospheric warming due to dominating non-monsoonal winter precipitation in this region. Fujita (2008) finds that HMA's glaciers are more affected by precipitation seasonality and concentration than by changes in annual precipitation. Where accumulation and warming happen at the same time (i.e. summer), rising temperatures increase both melt and the share of precipitation that falls as rain instead of snow. While temperatures are rising in entire HMA, the 10 glacier sensitivity study of Fujita and Nuimura (2011) suggests that temperature was not the limiting factor for glacier existence everywhere. In the extremely dry and cold TP and Kunlun Shan, with glaciers and in particular their accumulation areas at high elevations (Fig. 1), glacier growth due to increased precipitation is thus entirely plausible -despite a warming trend. This also stresses that the altitude of HMA glaciers (Fig. 1) is an important factor in their respective responses to temperature and precipitation changes (Sakai and Fujita, 2017), and thus in the here-observed glacier volume changes.

15
In light of continued climatic changes in the study region, ICESat only provides a short snapshot of ongoing glacier reactions.
This snapshot falls exactly into the decade where an increase in precipitation on the TP would cause the largest effects on glacier volume changes due to dynamic adjustment of the geometry in ablation areas as a delayed signal towards a new glacier equilibrium state Gilbert et al., 2018). dh-elevation gradients are largest, which could indicate that dynamical changes are happening also there: an overall thinning signal could be composed of increased melt at lower elevations, causing strongly negative dh, while the accumulation areas are growing due to increased precipitation/accumulation, causing strongly positive dh (Fig. 3d).

30
The negative trends on the eastern border of the TP agree with reported glacier mass loss in this area, although varying annually and in space (Kang et al., 2009;Yao et al., 2012). For this part of the southeastern TP, Mölg et al. (2014) found that the competition between the monsoon and large-scale westerly waves of the mid-latitude circulation in spring/early summer determines annual mass balance. The south-north transition of the jet stream across the TP in spring varies in timing and efficiency, and its re-intensification in summer on the northern edge of the TP is related to the onset of the summer monsoon (Schiemann et al., 2009). This interplay affects both precipitation and summer air temperature. All glaciers in the region are of summer accumulation type, except for East Nyainqêntanglha Shan and Hengduan Shan . The area where the atmospheric flow strength over the TP correlates strongly with summer temperatures  forms an arc-shaped band from the above mentioned mountain ranges along the northern slopes of the East Nyainqêntanglha Shan to

Glacier sensitivity to precipitation in the Himalaya
We find consistently less negative glacier volume changes on the first orographic ridge across the entire Himalayan Range.
Misclassifications of e.g. perennial snow patches with stable surface elevations classified as glaciers would cause a mixed glacier/land trend with a weaker surface lowering signal. To achieve this effect, the misclassification would have to be severe 25 (ca. half of the samples) and be present in both our manual classification and the RGI, as the pattern is visible with both glacier classifications. We carefully classify our samples manually to avoid precisely such mixed signals, thus we consider this bias unlikely. Another cause could be reduced melt due to insulation from debris cover. It has previously been shown that stagnant (debris-covered) tongues lose mass at a similar rate as clean ice glaciers Gardelle et al., 2012b;Pellicciotti et al., 2015;Ragettli et al., 2016). We thus assume that debris-cover is not the cause of the observed differences.

30
Locally varying sensitivity to precipitation might also explain the less negative mass balances on the first, and thus wettest, orographic ridge in the Himalaya. Precipitation from summer monsoon influx decreases sharply after large changes in relief (Bookhagen and Burbank, 2006). Maussion et al. (2014) find that precipitation regimes are strongly varying over short distances in the Himalaya, not least due to glacier orientation on the windward or lee side of the a mountain range. Wagnon et al. (2013) and Sherpa et al. (2016) mention the meteorologically exposed location of Mera glacier (4949-6420 m a.s.l.) in the Khumbu region, Nepal, as a possible explanation of its roughly stable mass balance since 2007 when in-situ measurements began. This stands in stark contrast to the considerable mass loss seen in Pokalde and Changri Nup glaciers only 30 km further north (the latter are also smaller and located at lower and thus warmer elevations, which likely contributes to these differences). In our ICESat zonation, these glaciers are located in units H1 (−0.12 ± 0.25 m a −1 ) and H2 (−0.50 ± 0.32 m a −1 ). Wagnon et al.

5
(2013) note that in the DEM differencing study of Gardelle et al. (2013), larger glaciers in the same range as Pokalde/Changri Nup also seem to experience more surface lowering than Mera glacier further south. The ELA sensitivity study of Fujita and Nuimura (2011) is too coarse to confirm orography-related spatial differences across previously been attributed to rapidly shrinking accumulation areas, seen in rising firn lines in Landsat images (Kääb et al., 2015, 20 area called Spiti Lahaul). Kääb et al. (2015) see the same pattern for the strongly negative glacier evolution in Nyainqêntanglha Shan/Hengduan, which has low-lying accumulation areas. Thus, once the accumulation area becomes too small or disappears entirely, also abundant or increasing precipitation cannot compensate for melt due to increased temperatures (Sakai and Fujita, 2017).

25
The zonation we present here is the result of an compromise between within-unit glacier similarity, representative sampling, and stable glacier surface change trends. In the Karakoram/Kunlun Shan area, this approach is clearly more appropriate than sample grouping into a regular grid. The latter results in large trend uncertainties (Fig. 2b) larger spatial units such as a regular grid, this local peculiarity is not visible. Whether such signal is representative for all glaciers in a unit or not would require complete geodetic analysis of all glaciers and also a longer time span.

Varied pattern in Tien Shan
Glacier evolution in the Tien Shan has shown a spatially diverse pattern already in the last decades of the 20th century (Narama et al., 2010;Farinotti et al., 2015). Together with contributions from northerly areas, the Westerlies are the source of precipita-the inner ranges are of the spring/summer accumulation type (Sorg et al., 2012). In the ( Kronenberg et al., 2016), and also DEM differencing/modelling studies in the area found similar values (Fujita and Nuimura, 2011;Shangguan et al., 2015;Barandun et al., 2018). Our zonation does not consider glacier aspects, which seem to play an important role in explaining glacier melt over this region (Farinotti et al., 2015). both areal extent and water surface levels. The latter work results in volume change time series of >1300 lakes, much more than available so far. In more detail, we conclude: -Only carefully delineated spatial units show local patterns of glacier change that are diluted or hidden if samples are gridded. On a larger scale, the pattern we find in this study agrees with previous regional estimates based on ICESatbut provides finer detail.

30
-The pattern of glacier changes is spatially very varied because glacier elevations, and their sensitivity to temperature and precipitation changes vary spatially (Sakai and Fujita, 2017;Kapnick et al., 2014). -Lake volume changes on the TP reflect a clear and comparably sudden increase of water availability from 1997 through ∼2010 for the northern and eastern TP, but only minor changes in the southwestern TP and Qilian Shan. The observed lake changes correspond to a precipitation-equivalent 6-7 mm a −1 for the northern TP and 25 mm a −1 for the eastern TP, 10 from decadal averages between the 1990s and 2000s. According to MERRA-2 reanalysis data, the change can exclusively be driven by increased summer precipitation of 34-100 mm decadal difference between the 1990s and 2000s. Decreasing potential evaporation from reduced wind speeds is also suggested to have in general contributed to lake growth (with uncertain timing though). Only in some areas increased influx from glacier mass loss should have contributed to lake growth as the zone of lake growths roughly coincides with the zone of positive glacier mass balances. The magnitude of lake volume change, glacier mass balance and precipitation changes agree within each other, considering also evaporation.
-Glaciers on the TP changed their geometry during 2003-2008. In the northeastern TP/western Kunlun Shan, upper glacier surface elevations were stable while tongues were growing. Further south/east, upper elevations were thickening while the tongues were thinning due to both increased accumulation and melt. The further southeast on the TP, the 20 stronger the glacier thinning rates. Glaciers in the Qilian Shan were only moderatly losing mass.
-Along the entire Himalayan Range, glaciers on the first orographic ridge were thinning less than those further back in a drier climate, likely due to abundant precipitation on the first ridge, which causes ELAs to be at lower elevations.
Precipitation and ELA gradients might be very steep in the outermost ridges of the Himalaya.
-Glaciers in the Tien Shan were thinning rather more than in other parts of HMA, in particular those in the transition 25 between the Tien Shan and Pamir mountains. There are exceptions to this general trend: glaciers in the central Borohoro range (at higher elevations) and on the northern slopes of the Tarim Basin were close to balance, possibly due to precipitation increase.
-From a methodological point of view, this work stretches the applicability and precision of ICESat-derived elevation changes in rough and glacierised terrain further than was the case for previous studies. We carefully examined the 30 influence of how spatial units are delineated to derive ICESat-based glacier change over HMA as well as a range of potential biases and error influences on the analyses.
While the glacier change pattern presented in this study is robust and well explained by glacier sensitivities to climate change, our unit boundaries might not match areas of consistent glacier changes everywhere, despite best efforts. Low ICESat sample density prohibits a further refinement in areas with sparse glacier coverage. Other remote sensing data with finer spatial reso-5 lution could improve the pattern -for example DEM differencing from ASTER stereo-imagery (Brun et al., 2017) and other spatially extensive data available for the last decades, or also ICESat-2, once this data becomes available.  (Schutz et al., 2005). Elevation samples are separated by ∼ 170 m along ground tracks/orbits but up to 75 km between orbit paths in HMA. The ground track pattern was not repeated exactly during each overpass, as the near-repeat orbit mode was not activated 20 at lower latitudes (Schutz et al., 2005). Rather, individual ground tracks lie as far as 2-3 km from the reference ground track in HMA. A direct comparison between ICESat elevations is thus difficult in the region. Instead, double-differencing techniques are applied, i.e. comparing ICESat elevations with a reference DEM to receive elevation differences and analysing their subsequent evolution over time Gardner et al., 2013;Neckel et al., 2014;Kääb et al., 2015;Ke et al., 2015a).
Here, we use GLAS/ICESat L2 Global Land Surface Altimetry HDF5 data (GLAH14, release 34) which is optimised for 25 land surfaces (Zwally et al., 2012). From comparison with reference DEMs, elevation uncertainty of GLAH14 data was found to be on the order of decimetres to metres in mountainous terrain in Norway (Treichler and Kääb, 2016). Elevation biases and inconsistencies throughout ICESat's lifetime are of centimetre to decimetre magnitude and thus negligible compared to uncertainties from the underlying terrain and biases in the reference DEM Treichler and Kääb, 2016).

30
The DEM from the Shuttle Radar Topography Mission (SRTM, Farr et al., 2007;Farr and Kobrick, 2000) is a consistent DEM in the HMA region. We used the C-band, non-void-filled SRTM DEM version at 3 arc-seconds resolution (SRTM3, corresponding to 92 m in y, and 66-82 m in x-direction at 45/28°N) which is accessible from the U.S. Geological Survey at https://dds.cr.usgs.gov/srtm. The SRTM DEM used here is a product of single-pass C-band SAR interferometry from images acquired on 11-22 February 2000 (Farr and Kobrick, 2000). SRTM DEM nominal vertical accuracy is of the order of metres (Rodriguez et al., 2006). Treichler and Kääb (2016) found spatially varying vertical offsets on the order of metres to decimetres in mountainous terrain in Norway. They attributed the vertical biases to the fact that the SRTM DEM is a composite from several individual images and overpasses, and likely processed in (unknown) spatial sub-units. Offsets caused by shifts of sub-units 5 were not removed by global DEM co-registration, but the bias/uncertainties caused by them are within the nominally stated accuracy. On glaciers, larger elevation uncertainties are to be expected due to penetration of the C-band signal into ice and, even more so, into snow. Also dry sedimentary soils may be subject to radar penetration. The penetration is estimated to be in the range of several metres for glaciers in HMA (Gardelle et al., 2012a;Kääb et al., 2012Kääb et al., , 2015. The vertical offsets from DEM shifts or penetration increase the uncertainty of surface elevation changes -possibly also for and Kääb, 2016. As an alternative elevation reference, we used the SRTM DEM at 1-arc-second resolution (SRTM1) from https://earthexplorer.usgs.gov. The 1-arc-second DEM has undergone fewer revisions than the 3-arc-second DEM, making the data not necessarily superior, and most data voids are filled in with other elevation data that have different time stamps.
We therefore excluded the data void areas contained in the 3-arc-second DEM version also in the SRTM1 DEM to ensure that 15 we only use original elevation data from February 2000. Further, we did not explore or use the recently published TanDEM-X global DEM as it was not available during our processing. It remains to be investigated how potential advantages of this DEM (larger coverage, less penetration than C-band) balance potential disadvantages (longer time difference to ICEat period, temporal inconsistency from stacking). Also due to temporal inconsistency and substantial voids, we did not use the ALOS PRISM World DEM (AW3D) or the WorldView satellite optical stereo HMA DEM.

A3 Precipitation data
As an estimate for regional and temporal precipitation patterns for for our study with respect to the lake volume changes investigated.
Further, we use in-situ data from the five westernmost meteorological stations on the TP and Kunlun Shan (Fig. 1), provided by the China Meteorological Science Data Sharing Service Network. The stations are located relatively close to the area with reported glacier mass gain (we are not aware of any meteorological measurements on the northeastern TP). The data includes daily measurements of precipitation, mean air temperature, and for the four stations on the southwestern TP also evaporation.

5
The Global Surface Water dataset (Pekel et al., 2016) is a classification of the entire Landsat archive into monthly and annual maps of surface water (https://global-surface-water.appspot.com). The data is available within Google Earth Engine (Gorelick et al., 2017). To map the changing extents of Tibetan lakes, we used the variable occurrence which provides the classes no data, no water, water (for both monthly/annual data), and seasonal water (for annual maps only). Coverage is nearly complete (>98%) starting from 2000 but considerably worse for some years of the 90s (Pekel et al., 2016, for our areas of interest: 10 20-75% no data pixels in 1990, 1991, 1995, 1997 and 1998). We follow the double-differencing method explained in further detail in Kääb et al. (2012) and Treichler and Kääb (2016), with special consideration of issues mentioned in the above introduction. ICESat data and individual SRTM DEM tiles were converted into the same geographical reference system, co-registered (Nuth and Kääb, 2011), and reference elevations for Per spatial unit, we estimate glacier surface elevation change by fitting a robust linear regression (which minimizes an iteratively weighted sum of squares) through individual dh . To test the sensitivity of biased dh at either end of the studied time period, we also compute a t-fit (Treichler and Kääb, 2016) and a non-parametric Theil-Sen linear regression (Theil, 1950;Sen, 1968). Both alternative robust fitting algorithms better fit our dh distribution and are commonly used for 30 datasets with large natural variability and measurement errors. Our final estimate per spatial unit corresponds to the average of the three trend methods.

B1 Zonation
ICESat data needs to be grouped into spatial units to fit surface elevation trends. The samples within each spatial unit need to reflect the glaciers in a representative way. This condition is easier to fulfil if the glaciers are similar to each other, including their 2003-2008 mass balances and their variations. We tested automated clustering methods from ICESat dh directly, but were not successful. We therefore preferred to delineate spatial units manually, considering topographic and climatic setting, 5 elevation, visual glacier appearance, and input from literature and discussions with experts. In particular, we paid attention to orographic barriers. The zonation we present here is thus the result of an iterative manual process of re-defining spatial units until they yielded statistically stable and robust glacier surface change estimates. While the procedure is based on carefully applied expert knowledge, we are fully aware that our zonation is eventually a subjective one and certainly open to discussion.
As a control approach, we applied the same gridding method as Kääb et al. (2012Kääb et al. ( , 2015 to the entire HMA.

B2 Glacier hypsometry
We compute the relationship between glacier dh and elevation (hereafter called dh-elevation gradient) by fitting a robust linear regression through individual glacier samples' dh vs. elevation. Greater radar penetration in the accumulation areas and more prominent melting of tongues steepen dh-elevation gradients (e.g. Vijay and Braun, 2016;Ragettli et al., 2016). It is therefore very important to ensure ICESat's elevation sampling is consistent through time and representative for glacier hypsometry (see 15 introduction). Our primary approach to improve sampling hypsometry is to enlarge spatial units, but in some areas this would have in turn led to considerably reduced glacier similarity within the unit. To account for these conflicting cases, we computed glaciers within each spatial analysis unit; (D) assigning weights to samples depending on their elevation so that they match the glacier hypsometry, i.e. analogue to C but without removing any samples.
Methods A and B are based on the method used in Kääb et al. (2012Kääb et al. ( , 2015. If ICESat consistently samples lower (or higher) elevations than the reference hypsometry, methods A and B will not correct for this -they only correct elevation-induced 25 bias relative to the mean sampled elevations of all campaigns. Methods C and D, however, adjust the hypsometry so that it should become representative for the glacier elevations in the unit. All four corrections are here applied to all units, and both for ice and land samples separately. Our 'standard method' for the final glacier elevation change estimates corresponds to the average of all hypsometry-correcting methods (A-D) and trend methods (robust, t-and Theil-Sen trends). Additionally, we also compute trends for only the upper/lower 50% glacier elevations as from RGI hypsometries (samples above/below the median 30 RGI glacier elevation in each unit). The latter analysis violates mass conservation and should thus not be interpreted in terms of mass balance, but rather, for instance, for changes in glacier elevation gradients (e.g. Brun et al., 2017;.

B3 Correction of vertical bias
Glacier elevation difference dh may be subject of vertical bias originating from local reference elevation bias or snow fall during the second part of the autumn 2008 campaign. Local vertical bias may result from inconsistent reference DEM age or production, tiling and tile/scene misregistration, or locally varying radar penetration (in case of the SRTM DEM). To remove this bias, we compute a per-glacier elevation correction cG corresponding to the median dh for each glacier, according to 5 the method described in Treichler and Kääb (2016). In that study, the correction successfully reconciled annual ICESat-based glacier elevation changes with mass balance time series from in-situ measurements. Also in the present study, cG-corrected dh (in combination with above hypsometry methods A-D) remove the effect of a varying spatial composition of elevation offsets. However, as the correction results in lower sample numbers and removes parts of the signal where some glaciers are only sampled in the beginning and some other glaciers only in the end of the ICESat acquisition period, the correction shows areal extent and lake levels from SRTM and ICESat data, respectively. Our method is in parts similar to the methods used by previous studies investigating lake volume changes on the TP from satellite data (Zhang et al., 2011;Kropáček et al., 2012;Song et al., 2013;Zhang et al., 2013;Song et al., 2015;Zhang et al., 2017, e.g., ) but the inclusion of a DEM for deriving shoreline elevations, and thus lake water levels, in addition to altimetry data enabled us to produce volume change time series for one order of magnitude more lakes than derived previously.
We apply our procedure to the 1364 endorheic lakes on the TP and in the Qaidam Basin ( Fig. 1) with a maximum lake extent of > 1 km 2 . We generated here our own lake database since we found that existing collections, such as the Global Lakes and Wetlands Database (Lehner and Döll, 2004), are lacking numerous lakes that likely only emerged during the last two decades.

5
Consulting satellite imagery like Landsat data, we manually adjusted our lake database to remove delta-like seasonal wetlands from water inflow on sloping terrain from the lake masks, we excluded non-endorheic lakes (visible outflow), and we excluded inundated areas affected by human interventions (e.g. for salt production) (133 wetlands not included in the above number).

C1 Uncertainties and filtering
To minimise the effect of uncertainties in or erroneous estimates for individual years, we analyse time series in a summarised way through regression over time and as decadal averages. Uncertainties associated with the lake data used include misclassification of water area in the Global Surface Water dataset (Pekel et al., 2016), lake surface elevation errors and local bias in the SRTM DEM, and bias in ICESat surface elevation measurements. For each lake and year, we compute the percentage of missing data (e.g. from cloud cover or classification voids), and years with < 95% of data coverage within the lake masks are excluded from further analyses. Lake time series that, after removing these years of insufficient coverage, do not contain 15 any data from the 1990s are excluded entirely. For ICESat-derived lake levels, only lakes with measurements from at least three laser footprints each from at least five years are considered. Despite the lake areas and surroundings being extremely flat, SRTM DEM cells indicate up to 10 m elevation differences between neighbouring cells in a seemingly random way, and the SRTM DEM turns out to be the dataset within our lake change analysis with the greatest uncertainties. For some lakes, SRTM DEM errors result even in negative area-lake-surface elevation relationships, i.e. lake shore elevations seemingly decrease 20 for expanding lake areas which is physically not plausible. We therefore excluded all lakes with either a negative area-lakeelevation relationship or where the 26-year linear trends for lake area and lake surface level do not have the same sign. This is done both for ICESat-and SRTM-derived lake level estimates. The overall error for a decadal average lake volume stage is estimated as the error of the mean, and for decadal differences propagated as the root of the sum of squares of the two errors (RSS).

C2 Endorheic basins
To estimate the lake water volume change in a way that can be related to glacier mass balances and precipitation changes (i.e. mm w.e. per m 2 ), we summarise and spatially distribute the water volume changes of all lakes within spatially confined basins. These basins are based on endorheic catchments, but because many catchments only contain a single lake and exact catchment areas are not well defined on the TP (e.g., in very flat areas), we manually controlled, adjusted and aggre-30 gated the endorheic catchments of the USGS HydroSHEDS dataset at 15 arcsec spatial resolution (Lehner and Döll, 2004, https://hydrosheds.cr.usgs.gov) to larger basins of comparable size and consisting of in average 5 catchments. We define the total lake area per catchment (and basins) as the sum of the 1990-2015 median lake area of all lakes within the spatial unit, also including the endorheic inundated areas confined by human infrastructure mentioned above, which are otherwise excluded from analyses. To compute total water volume change per catchment, we assume that lakes excluded from the analysis (see previous subsections) behaved the same way as the average of the lakes we have sufficient data for, and subsequently scale the 5 total volume change accordingly. For total water volume change from decadal averages, we compute the error as the sum of the errors of all individual lakes' volume change (see above), again scaled according to the share of total lake area we have sufficient data for. This conservative approach of adding errors (instead of root-sum-of-squares, RSS, for instance) includes as a worst case the full correlation of the behaviour of all contributing lakes. Representativeness of samples within spatial units is the key requirement for robust glacier thickening/thinning estimates.
However, we found that enlarging spatial units was not always the best remedy to ensure sample representativeness: In some areas this would have considerably reduced glacier similarity within the unit. Applying a regular grid can have this same effect.
Consequently, only carefully adapted zones can show local peculiarities that are otherwise diluted.
Especially for small units with few samples, careful consideration of how potentially biasing factors interplay is important.
Our use of four different methods to ensure correct hypsometry sampling makes the results very robust. The overall pattern is not affected by zonation, small changes in sample composition (RGI outlines), or reference DEM (here: SRTM1). Of all 5 corrections, the most essential requirement is therefore that the regional glacier hypsometry is sampled appropriately, also over time. Locally, however, the different methods and corrections can result in considerable differences between glacier thickening/thinning rates. Especially where ICESat data is used on a local scale or as input for modelling studies, we strongly recommend to carefully assess the difference between hypsometry corrections, the effect of our per-glacier correction cG, and the influence of snow cover, in order to ensure a representative estimate and appropriate uncertainty.

10
Our snow correction affects trends significantly. In southern Norway, the study region for which the correction was developed, it removed a positive land trend but did not affect the glacier trend (Treichler and Kääb, 2016