Recent glacier and lake changes in High Mountain Asia and their relation to precipitation changes

We present an updated, spatially resolved estimate of 2003–2008 glacier surface elevation changes for entire High Mountain Asia (HMA) from ICESat laser altimetry data. The results reveal a diverse pattern that is caused 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. Glaciers in the 5 Eastern Pamirs, Kunlun Shan and central TP were thickening by 0.1–0.7 m a−1, and 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. We attribute the glacier thickening signal to a step-increase ::::::: stepwise ::::::: increase in precipitation around :: ∼1997–2000 on the Tibetan Plateau (TP). The precipitation change 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 10 lake extents from Landsat data and shoreline elevations from ICESat and 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. Taking into account 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, 15 which can fully explain both lake growth, and glacier thickening (Kunlun Shan) or glacier geometry changes such as thinning tongues while upper glacier areas were thickening or stable (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, were thinning less than glaciers in the dryer climate of the inner ranges. Thinning rates in the Tien Shan vary spatially but are rather stronger 20 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 respectively (Yao et al., 2012;Bolch et al., 2012;Mukhopadhyay and Khan, 2014). On the very dry TP, glaciers occur only on the sparsely spread small mountain ranges.
In interplay with the Siberian High further north (Narama et al., 2010;Böhner, 2006), the Westerlies are the dominant source of moisture for the mountains surrounding the dry Tarim basin (at ca. 1000 m a.s.l.) -the Tien Shan to the north, and Kunlun Shan to the south (Ke et al., 2015;Yao et al., 2012). The mountain ranges at the eastern margins of the TP (Qilian Shan, 5 Hengduan Shan, Minya Gongga) are also influenced by the East Asian Monsoon (Yao et al., 2012;Li et al., 2015).
In both the monsoonal and westerly regimes, precipitation decreases northward . Depending on the regionally dominant source of moisture, glacier accumulation happens at different times of the year Maussion et al., 2014;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 .
For the HMA glaciers with predominant spring/summer accumulation, glacier accumulation and ablation happen at the same 5 time. Besides rising temperatures, recent studies suggest that climate change and altered circulation patterns affect radiation regimes and thus also glacier ablation in HMA through, e.g., changes in evapotranspiration or cloud cover (Forsythe et al., 2017;de Kok et al., 2018). The seasonal timing of snow accumulation on glaciers thus 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, which depends on continentality (Shi and Liu, 2000;Kuhle, 1990) and, on smaller spatial scale, on glacier 10 location on or behind a mountain range that acts as a primary orographic barrier or causes orographic convection. 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 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. 15 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, Kraaijenbrink et al. 2017. 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 Gardelle et al., 2012b;Pellicciotti et al., 2015). In this study, we distinguish thus not explicitly between debris-covered and debris-free glacier tongues. 20

Data and Methods
In this section we give a short overview of the data and methods used. Details can be found in the Appendix.

Data
For deriving repeat elevations on glaciers and lakes, we use data from the NASA Geoscience Laser Altimeter System (GLAS) aboard the Ice, Cloud and land Elevation Satellite (ICESat) that measured the Earth's surface elevations in two to three cam- 25 paigns per year from 2003 to 2009 (Zwally et al., 2012, GLAH14). The campaigns were flown in northern autumn (∼ October-November), winter (∼ March), and early summer (∼ June). (Appendix A1).
As reference DEM for our ICESat processing and to derive lake shoreline elevations we use the DEM from the Shuttle Radar Topography Mission (SRTM, Farr et al., 2007;Farr and Kobrick, 2000). 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 30 1-arc-second resolution (SRTM1). (Appendix A2) As an estimate for regional and temporal precipitation patterns for the years 1980-2015 we use data from the reanalysis products MERRA-2 (Gelaro et al., 2017) and ERA Interim (Dee et al., 2011). We use monthly summarised values of the variables total precipitation, snowfall and evaporation. The two chosen reanalysis products have previously been found to model precipitation comparatively well in our study area (Chen et al., 2019;Cuo and Zhang, 2017;Sun et al., 2018). (Appendix A3).
Further, we use in-situ data from the five westernmost meteorological stations on the TP and Kunlun Shan (Fig. 1), provided 5 by the China Meteorological Science Data Sharing Service Network. The data include daily measurements of precipitation, mean air temperature, and for the four stations on the southwestern TP also evaporation.
We extract repeat lake coverage from the Global Surface Water dataset (Pekel et al., 2016) that is a classification of the entire Landsat archive into monthly and annual maps of surface water. The data are available within Google Earth Engine. Spatial coverage is nearly complete (>98%) starting from 2000 but considerably worse for some years of the 1990s. (Appendix A4).

Methods for glacier volume change
We use surface elevation measurements from ICESat data points on glaciers and surrounding stable terrain and 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 (Appendix B). 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 15 the surface elevation has changed on average over the time period studied. We used only samples from ICESat's 2003ICESat's -2008 autumn campaigns, the season with least snow cover in entire HMA, to avoid bias from temporal variations in snow depths (see introduction). After filtering, 74'938 glacier samples and about ten times as many off-glacier samples remain. ICESat data need to be grouped into spatial units to receive surface elevation changes. The samples within each spatial unit need to reflect the glaciers in a representative way -which means that the spatial units need to be chosen such that they group glaciers 20 that are similar to each other in terms of climatic and topographical attributes, including their 2003-2008 mass balances and variations thereof. Previous studies have used regular grids, the RGI regions or their own arbitrary zonation. These do not necessarily fulfil the above requirements. We considered automated clustering methods to receive spatial units from ICESat dh directly, but were not successful. We therefore preferred to delineate spatial units manually, considering topographic and climatic setting, elevation, visual glacier appearance, and input from literature and discussions with experts (Appendix B1).

25
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 satisfied these criteria. After computing linear regressions on glacier dh, we split or merged some of the previously drawn units such that the final zonation 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 minimises an iteratively weighted sum of squares) and also compute a t-fit ) and a non-parametric Theil-Sen linear regression (Theil, 1950;Sen, 1968). Our 'standard method' for the final glacier elevation change estimates corresponds to the average of all hypsometry-correcting methods and linear regression methods. Additionally, we also compute elevation change for only the upper/lower 50% glacier elevations as from RGI hypsometries (samples above/below the median RGI glacier elevation of each individual glacier) for each spatial unit. The latter analysis violates mass conservation 5 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;Kääb et al., 2018). To allow comparison with other studies, we use RGI glacier areas to convert our surface elevation change rates to volume changes.
Glacier dh may be subject of vertical bias from elevation differences that are caused by other reasons than glacier surface elevation change, i.e. from bias in the local reference elevation (the SRTM DEM) or snow fall. We compute corrections for 10 these biases (Appendix B3). 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 the 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).

15
Analogue to Treichler and Kääb (2016), we estimate December 2008 snow bias from a linear regression of October/December 2008 off-glacier dh on elevation and time.

Methods for lake volume change
We derive the volume changes of lakes on the endorheic TP in order to relate glacier changes and precipitation changes on the Tibetan Plateau to each other. In particular, we want to investigate whether precipitation increases could be a reason for SRTM and ICESat data, respectively. Our method is in parts similar to the methods used by previous studies (e.g. Zhang et al., 2017) 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) than derived previously.
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).

5
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 aggregated from the endorheic catchments of the USGS HydroSHEDS dataset (Lehner and Döll, 2004). (Appendix C2).

Glacier thinning and thickening
10 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 orientation in order to group glaciers under similar temperature and precipitation regimes rather than across orographic barriers.
Surface elevation change for the new spatial units and the 2°× 2°grid in Fig. 2b  the other hand, reveal locally varying signals that are averaged out or not significant in the coarser grid zonation (e.g. units K1-K3 and KS1 in Fig. 2a). Our new zonation, surface elevation changes and corresponding glacier volume changes in Gt a −1 (using RGI glacier areas, Supplementary Information S2) are available as a data supplement.
In the Himalaya, the manual zone delineation shows a clear transition from moderately negative elevation change on the first, southern orographic ridge (−0.15 to −0.34 m a −1 , maximum trend error: 0.31 m a −1 ) compared to glaciers located further back The dh-elevation gradients are in some units very steep (Fig. 3d). This means that the surface elevation differences between ICESat and the SRTM DEM are very negative on glacier tongues but very small or even strongly positive in the upper accumulation areas. Steep dh-elevation gradients can result from altitudinal dependency of radar penetration or glacier geometry changes between SRTM and ICESat surface elevation acquisitions. The steeper the dh-elevation gradients are, the stronger is  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 15 correction to the error budget (Fig. 2c). The corrected glacier surface elevation change rates are included in the data supplement.
A discussion of the effect of this and other corrections and biasing influences is provided in Appendix D.

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 20 lake surface elevations (59% of the lake area). Extrapolated lake levels based on annual or campaign ICESat data (Appendix C) yield the same results, but ICESat-based lake level change is on average 1.55 times larger than SRTM-based values. Likely, the reason for this difference is the greater uncertainty of SRTM DEM elevations and pre-2000 SRTM lake levels (Appendix C1).
Multiplied with areal changes to receive volume changes, the relative difference is reduced to 1.09 times. Average 1990-2015 water level increase corresponds to 0.14 m a −1 (SRTM) and 0.18 m a −1 (ICESat) in lake-level change per year (  are unknown). Rather than growing steadily, most lakes seem to have undergone a phase of sudden and rapid growth starting between 1995 and 2000, and gradually slowing down until ∼2009, with rather stable conditions before and after this period.
(Note that lake time series are median-filtered due to data scarcity for the years 1995-1999. There is thus some uncertainty on the exact timing of the onset of lake growth.) Relative lake volume change was most sudden and rapid for the northeastern,  Endorheic lakes in the Qaidam basin/Qilian Shan 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).
5 Figure 4b shows the corresponding specific annual water volume change per endorheic catchment as the decadal difference between 1990-1999 and 2000-2009 average lake volumes (based on SRTM lake levels). The pattern of predominant water volume increase especially in the northern and eastern TP compares well to the results in Fig. 4a. Lake volume growth on the eastern TP is accentuated 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 Table 1. Water volume changes between decadal averages of the 1990s and the 2000s (dP and dV), and annual glacier mass balance of adjacent glacierised areas for 2003-2008 (last column). dV: total decadal lake water volume difference per basin region in mm m −2 , dP: annual precipitation difference in mm m −2 a −1 , station order in southwest TP: Shiquanhe, Gaize, Pulan, Nielaer. Glacier surface elevation changes are converted to mm w.eq. a −1 assuming a density of 850 kg m −3 .
Region dV SRTM dV ICESat dP MERRA-2 dP ERA Interim dP stations Glacier mass balance Northwest TP 62±14 70±15 34±11 −33 ± 11 16 ± 72 29 ± 10 to 31 ± 9 Northeast TP 60±12 54±9 85±13 −2 ± 22 13 ± 11 to 50 ± 21 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 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 5 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 also is plausible from Fig. 5.

Precipitation increase on the TP
A change in precipitation could explain both lake growth and glacier mass balance (if dominated by precipitation rather than temperature/melt). When subtracting the part that is lost through evaporation, the precipitation change should yield numbers that are directly scalable in relation to glacier mass balance and endorheic catchment water volume (when neglecting changes in subsurface water transport).

15
Annual precipitation sums on the TP from meteorological stations range from as little as 50-100 mm a −1 (Shiquanhe and Tashkurgan 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 southern- most stations close to the Himalaya (only ca. 50% precipitation in summer). On the data-sparse TP, both station data and reanalysis products may contain bias due to the stations not being representative for a larger area and the lack of observational forcing data for reanalysis products, respectively. We thus use the data in a summarised way and focus on relative changes rather than relying on absolute numbers to detect/confirm temporal changes and large-scale spatial patterns.
Of the five meteorological stations available, especially Shiquanhe and Pulan show little change in precipitation and pan 5 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 greatest relative change for the Gaize station (+42 mm per decade, a 25% increase) and Tashkurgan station (+16 mm or +22% 10 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 (Fig. 7), and also the spatial patterns of precipitation changes differ (Fig. 8). Figures 7a (ERA Interim) and 7b (MERRA-2) show regional averaged annual sums for total precipitation, evaporation, and the difference of the two, for grid points within the TP lake Winter precipitation did not change noticeably (−9 to −2 mm decadal change for the five TP regions, +8 mm for Qilian Shan). 5 Fig. 8b shows the spatial distribution of summer precipitation change (difference between decadal averages). Compared to the same map with ERA Interim data, MERRA-2 suggests a considerably stronger precipitation increase on the TP and increasing precipitation also in the Kunlun Shan area. For both reanalyses, the spatial patterns are the same for annual precipitation rather than summer precipitation only (not shown).
The two reanalysis products agree somewhat better when precipitation numbers are corrected with estimates of actual evap-10 oration to assess the total decadal increase in water availability. For MERRA-2, the decadal difference is then reduced to 6±11 mm (central TP) to 68±13 mm (northeastern TP). However, the evaporation-corrected increase is greater when looking at summer months only (31±7 mm to 77±11 mm per decade, compared to a decrease in water availability during winter months of −27±4 mm to −6±3 mm, not shown). Corresponding ERA Interim increase in annual water availability is 14±32 to 38±12 mm (summer: −19±11 mm in the Kunlun area to 38±13 mm, winter −16±6 to 5±4 mm). Both datasets suggest that 30-15 60% (MERRA-2) or 13-50% (ERA Interim) 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 the lower (TP) precipitation increase and, particularly, the decrease for the Kunlun Shan region (Fig. 8a) does not fit well with the results from our lake and glacier data.
Our above results on precipitation changes relate to decadal means in order to enable systematic comparison to other data. It 5 is however important to note that these results vary if other time periods are chosen for aggregation. Kääb et al. (2018), 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 periods.

Discussion
The 2003-2008 ICESat surface elevation changes paint a spatially diverse picture of glacier changes in HMA. The general 10 pattern -glacier volume gain in the Kunlun Shan and the inner TP and glacier volume loss elsewhere -appears robust, no matter whether we aggregate the samples in a regular grid or manually delineated units. The more distinct spatial pattern agrees with the ICESat studies of Kääb et al. (2015Kääb et al. ( , 2012, the ASTER-based geodetic mass balances of Brun et al. (2017) and with the overall picture drawn by the previous regional studies of Neckel et al. (2014), Gardner et al. (2013) and Farinotti et al.
(2015) based on data from ICESat, GRACE and modelling. The pattern found is also robust against small changes in reference 15 elevations (such as from using the 1 arc-second SRTM DEM) or sample composition, and can also be reproduced using the most recent RGI glacier outlines -which have clearly become much more accurate since the study of Gardner et al. (2013).
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 20 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. A method discussion, in particular on biasing influences on ICESat glacier surface elevation change rates, is provided in Appendix D.

Coincident lake growth and glacier thickening 25
The regions with glacier thickening, or thickening of upper glacier areas (Fig. 3a), spatially match the areas with growing Landsat data gaps (1996,1997, see App. A4). The recent growth of TP's lakes is established by numerous recent studies (e.g., Zhang et al., 2011Zhang et al., , 2013Song et al., 2015;Zhang et al., 2017). In this study, lake volume changes on the TP serve as proxies for precipitation changes, but they may also help resolving satellite gravimetric signals to compute glacier mass changes (see introduction). The fact that glacier volumes are increasing in regions where also lake volumes increase, and the fact that lake volumes are also increasing in little or not glacierised 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).

5
Though, glacier mass loss can certainly play an additional role for lakes with declining glaciers in their catchment. This is in line with, and extends geographically, water balance modelling by Lei et al. (2013) for six selected lakes in our East TP zone (Fig. 4b) 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 (Song et al., 2015). Evaporation may also have decreased due to increased humidity from higher precipitation amounts. For 1981-2013, roughly stable (MERRA-2) evaporation for the two regions in the north. It is noteworthy that correcting precipitation data with evaporation allows to somewhat reconcile the two reanalysis datasets: Also ERA Interim shows an increase in so-computed net water availability, although it is smaller than for MERRA-2. Lei et al. (2013) suggest that groundwater exchange between different basins has very limited influence on the water balance 20 of each lake due to the impermeability of surrounding permafrost. Such groundwater exchange does not affect the basinwide water volume changes of this study, but thawing permafrost could be another potential source of water. An increase of active layer depth also causes an increase in groundwater storage capacity in ice-free ground an may change the amount of precipitation or water from snow melt that is retained or released (S. Westermann, pers. comm.). However, we are not aware of studies that quantify the amount of water available from these processes. Modelling studies (Ran et al., 2018;Zou et al., 2017) 25 find continuous permafrost in the northern part of the TP (our regions NW, NE, C and most of E) and discontinuous permafrost including larger areas of non-frozen ground in the southern/eastern parts of the TP (our regions SW and most of QQ). Recent and ongoing temperature rise led to an increase in the active layer and degrading permafrost that seems to have been greatest during the 60s and 00s and in the southern and eastern parts of the TP (Ran et al., 2018), where we find little lake change / lake growth (SW) and strong lake growth (E), respectively.

Precipitation increase on the TP and glacier sensitivity to these changes
In particular the MERRA-2 reanalysis, and to a lesser degree also ERA Interim and station data, suggest precipitation on the TP has increased around 1995-2000. The spatial patterns of decadal precipitation increases and the glacier growth on the TP and in the Kunlun Shan suggest a causal relationship. Increased precipitation in the region has been noted before: Yao et al. (2012) attributed a pattern of precipitation/glacier changes to a strengthening of the Westerlies while the Indian monsoon is weakening.
A rise in extreme precipitation events at stations in the study region was attributed to a weakening East Asian monsoon (Sun and Zhang, 2017). Fujita and Nuimura (2011) and Sakai and Fujita (2017)  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 glacier sensitivity study of Fujita and Nuimura (2011) suggests that temperature was not the limiting factor for glacier existence 10 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 elevation 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.

Glacier geometry changes on the TP
In The rate of warming on the TP is greatest for the elevations where glaciers have their ablation areas (Yao et al., 2012;Ran et al., 2018). In the southeastern part of the TP, dh-elevation gradients are largest (darker units in Fig. 3d), which could indicate that dynamical changes are happening also there: an overall thinning signal could be composed of increased melt 30 at lower elevations, causing strongly negative dh, while the accumulation areas are thickening or stable due to increased precipitation/accumulation, causing stable surface elevations or positive dh. This interpretation is supported by the gradual transition visible in Fig. 3c: in East Kunlun Shan and central TP, we see a thickening of accumulation areas and no change on the tongues (area marked "G"), and further east/south accumulation areas experienced little change but tongues were thinning (marked "L/A").
Dynamic glacier geometry adjustments might also be reflected in glacier flow. Dehecq et al. (2019) found that, for the 2000-2016 period, the flow speed of HMA glacier tongues decreased everywhere but in the Kunlun Shan and Karakoram and only slightly decreased on the TP. While the different time periods and spatial aggregation don't allow a more detailed comparison, 5 their results confirm that these regions were not or less affected from rapid glacier mass loss with thinning and increasingly inactive tongues.

Glacier thinning on the Eastern Tibetan Plateau
The negative elevation change rates on the eastern border of the TP agree with reported glacier mass loss in this area, although Glacier lies, and even less negative values further west (−0.14 ± 0.10 m w.e. a −1 ), in line with Tian et al. (2014). Towards east, glaciers become smaller and elevations lower, and the influence of the East Asian Monsoon becomes stronger.

Glacier mass balance and precipitation in the Himalayas
We find consistently less severe glacier thinning 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 (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. 5 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.
A potential explanation for the less negative mass balances on the first, and thus wettest, orographic ridge in the Himalaya 10 is a locally lower sensitivity of glacier mass balances to precipitation (and changes thereof). 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 The ELA sensitivity study of Fujita and Nuimura (2011) is too coarse to confirm orography-related spatial differences across the Himalaya, but along the mountain ridge their findings correlate well with both Yao et al. (2012) and our pattern of glacier changes in the inner Himalayan ranges (see also Sakai and Fujita, 2017). In particular the stable glacier elevations in our unit 30 HK1 -between areas of glacier loss in the Hindu Kush and the particularly negative western Himalaya (units H4-H6)are backed up by their modelled stable ELAs. According to MERRA-2 data (but not ERA Interim), the area experienced an increase in summer precipitation between the 90s and 00s (Figs 8a, 8b). The particularly negative surface elevation change in the western Himalaya has previously been attributed to rapidly shrinking accumulation areas, seen in rising firn lines in Landsat images (Kääb et al., 2015, 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).

5
The zonation we present here is the result of an compromise between within-unit glacier similarity, representative sampling, and stable glacier surface change rates. In the Karakoram/Kunlun Shan area, this approach is clearly more appropriate than sample grouping into a regular grid. The latter results in large uncertainties in the glacier elevation change signal (Fig. 2b), since grid cells include both the thinning signal south of the central Karakoram and thickening signal in the Kunlun Shan.
In the Karakoram, we see indications of both surging glaciers and glaciers recovering from a surge. In most units, such as

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 precipitation for the entire region (Bothe et al., 2012), but there are different climatic sub-regions: glaciers in the Western Tien Shan (and Pamir Alai) receive precipitation mainly in winter, the northern and northeastern ranges both in winter and summer, whereas 25 the inner ranges are of the spring/summer accumulation type (Sorg et al., 2012). In the (north)western Tien Shan, our zonation does not consider this transition from winter-only to summer/winter precipitation due to too low sample numbers for a finer zonation in this area. Narama et al. (2010) suggest that glaciers of the outer ranges -which receive more precipitation -are melting faster since they have a higher mass turnover and their tongues are at lower elevations. They see such a pattern in 2000-2007 glacier 30 shrinkage, which was more pronounced in the Western/northern Tien Shan than in interior areas such as the southeastern Fergana Range or At-Bashy Range at the transition to the Pamir. Our thinning rates do not confirm this -precisely in this latter area (unit P3), we find the most negative glacier surface elevation changes in the entire region (converted to mass change: −1.04 ± 0.23 m w.e.). The modelling study of Farinotti et al. (2015) suggests spatially highly varying glacier reactions in the last few decades in that area (their coarser zonation in the Central Tien Shan does not allow direct numerical comparison with our results).
ICESat suggests moderate thinning for the north-eastern Borohoro range, in particular the central part at higher elevations (TS1, converted: −0.09 ± 0.18 m w.e. a −1 , upper 50% glacier elevations thickening in Fig. 3a). Farinotti et al. (2015) found that the central parts of the range receive 50% more summer precipitation compared to the rest of the range, and modelled −0.17 ± 0.24 m w.e. a −1 for 2003-2009 for a slightly larger area than our most central unit.
In the inner Tien Shan, our elevation change rates vary on a small spatial scale. Reconstructed annual mass balances (Kenzhebaev et al., 2017;Kronenberg et al., 2016) and DEM differencing/modelling studies in the area (Fujita and Nuimura, 2011;Shangguan et al., 2015;Barandun et al., 2018) match the range of our thinning signal. Our zonation does not consider glacier 10 aspects, which seem to play an important role in explaining glacier melt over this region (Farinotti et al., 2015). gridded. On a larger scale, the pattern we find in this study agrees with previous regional estimates based on ICESat -but provides finer detail. The new zonation and improved bias control in 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.
-The pattern of glacier changes is spatially varied because of differences in the glaciers' elevations and sensitivity to climate changes (Sakai and Fujita, 2017;Kapnick et al., 2014). Together with glacier elevations, precipitation distribution and changes are able to explain large parts of the spatial variability of the glacier change pattern observed for 2003-2008.

5
-An almost step-like precipitation increase on the TP, Kunlun Shan and possibly also the Tarim Basin between 1995 and 2000 is clearly visible from changes in lake water volume as well as MERRA-2 reanalysis data. The precipitation increase is able to fully explain 2003-2008 glacier thickening in an area centred over the Kunlun Shan. The boundary between positive and negative glacier changes is rather sharp in the Kunlun Shan and lies north of the main Karakoram range. It is more gradual on the TP. Also glaciers on the northern slopes of the Tarim Basin were close to balance.

10
-Lake volume changes on the TP reflect a clear and comparably sudden increase of water availability from ca. 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, from decadal averages between the 1990s and 2000s. MERRA-2 reanalysis data suggests the change is exclusively be driven by increased summer precipitation of 34-100 mm decadal difference between the 1990s and 15 2000s. ERA Interim reanalysis data suggests a smaller precipitation increase for a smaller spatial area that does not explain lake growth and glacier thickening equally well.
-The magnitude of lake volume change, glacier mass balance and precipitation changes agree within each other, considering also evaporation. Increased influx from glacier mass loss may in some areas have contributed to lake growth but cannot explain it, as the zone of lake growths roughly coincides with the zone of positive glacier mass balances or 20 dynamical glacier geometry change.
-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 stronger the glacier thinning rates. Glaciers in the Qilian Shan were only moderatly losing mass.

25
-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 equilibrium line altitudes (ELAs) to be at lower elevations. Precipitation and ELA gradients might be very steep in the outermost ridges of the Himalaya.
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 sam-30 ple density prohibits a further refinement in areas with sparse glacier coverage. Other remote sensing data with finer spatial resolution 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. Combinations of remote sensing products for precipitation, snow and atmospheric parameters as well as improved reanalysis data could help to determine precipitation numbers with more certainty in Asia's water tower.
Code and data availability. ICESat data are freely available from NSIDC and NASA, the SRTM DEM and Landsat data from USGS, the MERRA-2 reanalysis data from NASA Goddard Earth Sciences Data and Information Services Center, ERA Interim from the European 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 were not removed by global DEM co-registration, but the bias/uncertainties caused by them are within the nominally stated 5 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 ICESat-based studies, if the spatial pattern of SRTM DEM offsets interferes with ICESat's spatial sampling pattern (Treichler 10 and . 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 we only use original elevation data from February 2000.

15
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 ICESat 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 the years 1980-2015 we use data from the Modern-Era products have been found to model precipitation and snowfall comparatively well (Reichle et al., 2017a, b). The High Asia Reanalysis (HAR, Maussion et al., 2014), 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.
The meteorological stations included in this study were chosen because they are closest to the area with reported glacier mass gain. We are not aware of any meteorological measurements on the northwestern TP.

A4 Global Surface Water Dataset
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 5 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). Pre-2000 coverage is poor for years with little Landsat data, for our areas of interest: 20-75% no data pixels in 199020-75% no data pixels in , 199120-75% no data pixels in , 199520-75% no data pixels in , 199720-75% no data pixels in and 199820-75% no data pixels in (Pekel et al., 2016. Appendix B: Methods for glacier volume change 10 We follow the double-differencing method explained in Kääb et al. (2012) and Treichler and Kääb (2016). 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 ICESat footprint centres retrieved by bilinear interpolation. 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 the surface elevation has changed on average over the time period studied. To test the sensitivity of biased dh at either end of the studied time period, we do not only compute a robust linear regression, which is commonly used for ICESat glacier applications , but also a t-fit ) and a non-parametric Theil-Sen linear regression (Theil, 1950;Sen, 1968). Both alternative robust fitting algorithms better fit our dh 25 distribution and are commonly used for datasets with large natural variability and measurement errors.
We find little difference between robust and t-fits, and slightly larger (but no systematic) differences when using Theil-Sen linear regression. The trend slopes from the three methods agree on average within 0.1 m/a and differences always lie well within trend error estimates. Our final estimate per spatial unit thus corresponds to the average of the three trend methods.

B1 Zonation
As seen in Kääb et al. (2015), grouping of ICESat samples into a regular grid without a-priori knowledge results in a blurring of local glacier change signals. Since such local signals consist of a specific dh magnitude and evolution over time which should be governed by climatic or topographic drivers, we tried to derive a more realistic spatial division from the ICESat samples directly, using glacier statistics, dh and iterative clustering. This approach was not successful: the number of (semi-5 quantitative) statistical parameters turned out to be too large and dh vary too much spatially, not least due to bias. We thus carefully delineated spatial units manually. Zones were drawn by hand to avoid splitting any glacier between several zones.
In particular, we paid special attention to orographic barriers. Rather than roundish zones across the entire Himalayan range, we chose elongated zones around mountain ridges. Size, length and width of spatial units (i.e. how many parallel ridges) were largely determined by ICESat sample numbers and the condition of representativeness. For example, we included both the 10 windward and leeward side of a Himalayan range as there are very few glacier facing south (i.e. windward), and we suspected that leeward accumulation areas close to a mountain peak might still receive more precipitation from turbulences than the dry, leeward valley bottoms (Immerzeel et al., 2014). We are very aware that our zonation is a subjective one and open to discussion.
In some parts, other operators will likely come up with modified zones. However, our zonation is based on carefully-applied expert knowledge, and we are convinced it displays the 2003-2008 HMA glacier elevation changes with a spatial resolution 15 and precision that reflects the optimum that is feasible from ICESat over such a mountainous and heterogeneous region.

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). Representative 20 elevation sampling through time and in relation to local glacier hypsometry is thus very important. Our primary approach to improve sampling hypsometry is to enlarge spatial units, but in some areas this would have led to considerably reduced glacier similarity within the unit. To account for these conflicting cases, we computed four different corrections and compared the such- (C) filtering of the samples of each ICESat campaign to match the hypsometry of the 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.

30
All four corrections are here applied to all units, and both for glacier and off-glacier samples separately. 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 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. For 18 units, the difference in derived surface elevation change between the 'standard method' (average of all methods A-D) and only applying the latter methods (average of methods C and D) exceeds 0.05 m a −1 , and at the same time, average glacier elevations sampled in these units is also > 50 higher or lower than average glacier elevations for this unit (SRTM elevations within RGI glacier outlines). For these units with systematic elevation 5 missampling, we used the average of methods C and D only. To 5 of the affected units we also applied the cG correction (see below).

B3 Correction of vertical bias
To remove local systematic elevation bias, we compute a per-glacier elevation correction cG corresponding to the median dh for each glacier (i.e. subtracting the median dh for each glacier from each corresponding dh). In the study of Treichler and  . We estimate the influence of this according to the method of Treichler and Kääb (2016).
Appendix C: Methods for lake volume change 25 We compute annual water volume change of the Tibetan lakes by multiplying annual lake areas with water level changes 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 by exporting bitmaps of annual water occurrence over the entire TP, using the web API of Google Earth Engine. The data is exported at a resolution of 50 m × 38-44 m in lat/lon (corresponding to 0.00045 degrees).
Subsequently, we retrieve the corresponding lake surface elevations in two ways: a) from SRTM DEM elevations of the lake 30 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 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.
We extract the relationship both for annual time series and individual ICESat campaigns (2-3 campaigns each year, using the monthly water classifications). 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 5 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 10 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. 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 (in total 133 wetlands not included in the 15 above number). For spatial aggregations, computation of relative numbers per lake and for plotting, we use the median lake areas from the 1990-2015 annual lake extents.

C1 Uncertainties and filtering
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 20 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 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. Data from the 90s have higher uncertainties in extracted/extrapolated lake levels due to a) the implicitly assumed 25 bathymetric profile using area-lake level scaling for years without ICESat data; and b) because the SRTM DEM was acquired in February 2000: While lake areas vary seasonally and we use annual maximum areas, the effect of extracting SRTM lake elevations for lake areas smaller than during SRTM data acquisition is that some pre-2000 SRTM lake levels may be too high, resulting in too small dV. 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 30 within our lake change analysis with the greatest uncertainties. Potential explanations for the DEM elevation uncertainties are penetration of C-band radar into sandy ground and unknown processing steps during DEM production to mask/interpolate water-covered areas without radar backscatter. For some lakes, SRTM DEM errors result even in negative area-lake-surface elevation relationships, i.e. lake shore elevations seemingly decrease for expanding lake areas which is physically not plausible.
We therefore excluded all lakes with either a negative area-lake-elevation 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 standard 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 5
We summarise and spatially distribute the water volume changes based on endorheic catchments of the USGS HydroSHEDS dataset at 15 arcsec spatial resolution (Lehner and Döll, 2004, https://hydrosheds.cr.usgs.gov) . However, many catchments only contain a single lake and exact catchment areas are not well defined on the TP (e.g., in very flat areas, Lehner et al., 2008) and the spatial resolution of the HydroSHEDS dataset is in parts too coarse to correctly attribute the lakes of our lake dataset to the correct catchment. Therefore, we manually controlled and adjusted the endorheic catchment borders using the finer 10 topography of the SRTM DEM at 3 arcsec resolution as well as Landsat imagery to detect surface water exchange between lakes/catchments, and aggregated the catchments to larger basins of comparable size, 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 15 the analysis (see previous subsections) behaved the same way as the average of the lakes we have sufficient data for, and subsequently scale the 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. 20 Appendix D: Discussion of biasing influences on ICESat glacier surface elevation change 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 the same effect.
Consequently, only carefully adapted zones can show local peculiarities that are otherwise diluted.

25
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 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 30 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.
Our snow correction affects trends significantly. In southern Norway, the study region for which the correction was developed, it removed a positive off-glacier trend but did not affect the glacier trend . Our results in HMA show that trend fitting methods are surprisingly sensitive to a lowering of the last (half) campaign, no matter which trend fitting 5 algorithm is used, and for both off-glacier and glacier dh. In contrast, if the same correction is applied to a campaign between 2004 and 2007, trends only change marginally. The exercise shows that November/December 2008 snow fall has the potential to erroneously decrease ICESat-derived glacier thinning rates, in particular in Tien Shan, Pamir, Hindu Kush, Nyainqêntanglha Shan/Hengduan Shan and maybe also the outer Himalayan ridges (Supplementary Fig. S1). We therefore recommend to assess the bias potential of December 2008/October snow fall for ICESat studies on a smaller spatial scale. Also, we advise not to rely on ICESat's March campaigns for glacier studies wherever snow is falling in winter in the northern hemisphere.
ICESat elevations have previously been used to estimate SRTM penetration Shangguan et al., 2015). On glaciers where no ICESat data is available, dh-elevation gradients of larger spatial units -such as in this studycould improve the estimated elevation dependency of penetration.