Hot Spots of Glacier Mass Balance Variability in Central Asia

Abstract The Tien Shan and Pamir mountains host over 28,000 glaciers providing essential water resources for increasing water demand in Central Asia. A disequilibrium between glaciers and climate affects meltwater release to Central Asian rivers, challenging the region's water availability. Previous research has neglected temporal variability. We present glacier mass balance estimates based on transient snowline and geodetic surveys with unprecedented spatiotemporal resolution from 1999/00 to 2017/18. Our results reveal spatiotemporal heterogeneity characterized by two mass balance clusters: (a) positive, low variability, and (b) negative, high variability. This translates into variable glacial meltwater release (≈1–16%) of annual river runoff for two watersheds. Our study reveals more complex climate forcing‐runoff responses and importance of glacial meltwater variability for the region than suggested previously.

and traced the albedo change with elevation. The greatest slope on the albedo-elevation profile defined the critical albedo (α crit ), implying the limit between snow and ice / firn, thus the transient snowline. The critical albedo was considered to be a site-and scenespecific albedo threshold (α crit ) and was used as reference value to re-evaluate the pixels within the range of critical albedo values of the specific scene. In order to reclassify extreme outliers, all grid cells were re-examined regarding their relative position compared to the beforehand defined snowline. We used a probability-based approach: an increasing positive / negative vertical distance from the transient snowline resulted in a decrease in likelihood of the cell to associate to the class on opposite site of the transient snowline.
Pixels were reclassified accordingly.
We derived the snow-covered area fraction (SCAF) i.e. the ratio of the area above the current snowline to the total area of the glacier, to include the glacier-specific hypsometry. Glacier outlines were taken from Randolph Glacier Inventory version 6 (RGI 6.0, RGI Consortium (2017)). Additionally, to avoid misclassification of completely snowcovered glacier surfaces, we assumed that the median albedo within the lower quarter of the glacier only exceeded the value of >0.35 in the case of an entirely snow-covered glacier. In a last step, the transient snowlines used for model calibration were filtered to avoid misclassification due to sporadic fresh snowfall based on the temporal evolution of the snowline throughout the ablation season. The minimum observed snow covered area fraction (SCAF) at the end of the summer defined the end of the ablation season, and all following SCAFs were omitted.
May 6, 2021, 12:55pm : X -5 The separation of snow and ice primarily depends on the albedo of the surface classes.
It is possible to account for fresh snowfall or changes in other surface characteristics with a physical interpretation of this value. However, the range of albedo for snow and ice / firn can overlap considerably, and a straightforward classification is not always possible (Naegeli et al., 2019). Nonetheless, the transition of the albedo between ice and snow is spatially characterised by a distinct change (Hall et al., 1987), and the greatest slope of the albedo-elevation profile approximates the limits between snow and bare-ice surfaces well (Naegeli et al., 2019). On average, 5-10 SCAFs are available per ablation season and glacier.
Comparison to manually delineated snowlines for Abramov, Golubin and Glacier No. 354 showed satisfying agreement within an RMSE of less than 10% for all mutually observed SCAFs. However, the validation focused on clear-sky images. The applied cloud mask using SAM does not filter very thin optical clouds and hence can interfere with automated surface classification.
Misinterpretation of the snowline on the satellite images can affect the model performance and a thorough filtering of snowlines is recommended. Barandun et al. (2018) found an uncertainty of ≈0.10 m w.e. yr −1 related to a systematic over-and underestimation of the mapped snowlines and showed that the model sensitivity to single snowline delineation is not strong if enough and well-distributed information on the transient snowline position is available. This highlights the benefit of the used method to not solely depend on the end-of-summer snowline as a proxy for the surface mass balance. Including all available X -6 : transient snowline observations hence helps to guarantee a robust calibration procedure.

Geodetic volume change
We produced ASTER DEMs using the MMASTER processing chain (Girod et al., 2017).
To remove biases due to sensor motion in the ASTER DEMs, we used the SRTM DEM as a reference, as it provides a consistent DEM for the entire region. We then masked the final DEMs using the correlation score provided by the MMASTER processing, using a threshold for the correlation score of 70%.
To filter the DEM differencing pairs, we employ a filtering procedure, modified after Pieczonka and Bolch (2015), to remove remaining outliers on the glacier area. The introduced filter used the overall standard deviation (σ) of the glacier elevation difference, weighted by an elevation-dependent coefficient, allowing both positive and negative elevation changes in the ablation and accumulation area. The weighted coefficient for each pixel was determined using a sigmoid function. After a normalisation of all pixels on the glacier surface to its elevation range (Eq. 1), a weighting coefficient c was calculated using the normalised glacier elevation w (Eq. 2). This coefficient was then multiplied with the standard deviation of the glacier elevation change σ dh to indicate the expected maximum change of each pixel on the glacier surface δh max .
where z x,y is the pixel on the glacier surface, z min and z max the minimum and maximum glacier elevation, respectively.
Because of low correlation due to clouds, shadow, or bright snow, the DEMs (and therefore the elevation change maps) had significant voids over the glaciers, which prevented a straightforward summation of the elevation differences to calculate volume changes.
Therefore, we used the local mean hypsometric approach  to calculate volume changes. First, we estimated the area-altitude distribution for each glacier based on the SRTM DEM, using either 50 m elevation bins or bins corresponding to 10% of the glacier elevation range, whichever was smaller. For each DEM differencing scene, we then took the mean elevation difference per elevation bin, iteratively removing remaining outliers that were more than three standard deviations away from the mean of the bin, which yielded an estimation of the elevation difference as a function of elevation. We then removed all values where the bin was more than 40% void (i.e., no data available). If the elevation curve still covered at least 75% of the glacier area based on the area-altitude distribution, the curve was filled using linear interpolation, setting elevation difference values above and below the glacier elevation range to zero. Elevation difference maps where less than 50% of the glacier area was sampled or the root mean square of the elevation difference on stable ground was more than 10 m were removed. Finally, we calculated The difference in volume changes ∆V for each glacier based on ASTER scenes were converted into geodetic mass balances ∆M geod , using a density ρ ∆V of 850 kg m −3 (Huss, 2013): where, A is the glacier area and ∆t is the time in years between the corresponding image pairs.

Homogenisation of geodetic mass balances
A large amount of overlapping ASTER scenes and HMA DEM for the studied glaciers led to a range of different geodetic mass balances with different time stamps per glacier (S1). In order to provide a suitable and robust result for a second-order calibration all geodetic mass changes for each individual glacier i were homogenised to represent a fixed reference period t ref (1999/00 to 2017/18). To homogenise geodetic mass balances covering arbitrary periods t to the reference period we followed Zemp et al. (2019). First, we calculated the mean annual deviation β t between each geodetic estimate M geod and : X -9 the selected glaciological time series M glac over a common time period of N years between t 0 and t 1 : In a second step, we calculated the mean annual mass change ∆M ref for t ref for each individual glacier by adding β t to the mean annual glaciological mass balance for the reference period M glac,t ref : For the homogenisation step, the series of Tuyuksu was used for Dzhungarsky Alatau,

Western / Northern Tien Shan, Pamir-Alay and Western Pamir and the one of Urumqi for
Eastern and Central Tien Shan and Eastern Pamir. All ∆M t ref ,i were weighted according to their specific uncertainty and the median of all the weighted estimates was interpreted as the reference geodetic mass balance of the corresponding glacier (S1). This value was then used for a second-order model calibration. For each glacier, the calculated uncertainties of the mass changes obtained for the different periods are summarised through their arithmetic mean. This value is used to represent the uncertainty of the geodetic mass balance for the reference period.

Comparing homogenised geodetic mass balances with other studies
A glacier-by-glacier comparison of the mass balance with results provided in Brun, Berthier, Wagnon, Kääb, and Treichler (2017) agreed within the expected uncertainties with a mean difference of −0.15 m w.e. yr −1 , ours being more negative, and RMSE X -10 : of 0.27 m w.e. yr −1 . Good agreement was obtained for the region-wide mean mass balance rate for the Tien Shan (<−0.05 m w.e. yr −1 ). Region-wide differences were larger for the Pamir (<−0.30 m w.e. yr −1 ). We compared the homogenised geodetic mass balances from 1999/00 to 2017/18 with the geodetic surveys published by Shean et al. (2020) and Hugonnet et al. (2021). We found good agreement for the Tien Shan ( Fig. S2 and Table S1). However, geodetic mass balance estimates for the Pamir presented in this study were somewhat more negative than in Shean et al. (2020) (Brun et al., 2017;Shean et al., 2020;Hugonnet et al., 2021) are mostly based on the same ASTER datasets. Independent mass change assessments for the Pamir are partly inconsistent however (Gardelle et al., 2013;Brun et al., 2017;Shean et al., 2020;Gardner et al., 2013;Kääb et al., 2015;Farinotti et al., 2015;Barandun et al., 2020), with our own estimate falling somewhere in between. In a first step, we constrained the model parameters C prec and DDF snow for each glacier and year with SCAFs derived from the transient snowline observations. In Barandun et al. (2018), an initial range for DDF snow and C prec was narrowed down until an optimal parameter combination was obtained. In order to limit computation time, the calibration procedure was adjusted for a regional application (Fig. S6).

First-order model calibration
C prec constrained the modelled accumulation of a specific year and glacier. It was adjusted to minimise the root mean square error (RMSE) between the modelled cumulative melt M elt modelled at the position of the observed snowlines, and the total amount of accumulated winter snow Accu modelled that melted from the onset of the ablation season until the snowline observation date (Fig. S7, Barandun et al. (2018)). DDF snow was calibrated to best represent all SCAF observations of one ablation season by reducing the RMSE between modelled and observed SCAF RM SE SCAF (Barandun et al., 2018). C prec and DDF snow were calibrated annually and for each glacier separately to correctly represent the winter snow accumulation and the melt rates for each year of the observation period.
More details are given in (Barandun et al., 2018).
We do not explicitly describe internal and basal balance and sublimation but account for them within the expected uncertainties by integrating geodetic mass balances into the calibration procedure. As a control, the mountain range mass balance of both methods (transient snowline constrained mass balance reconstruction / geodetic method) was compared and no systematic offset (mean absolute error = +0.09 m w.e. yr −1 ) was found : X -13 (correlation coefficient = 0.71; RMSE = 0.28 m w.e. yr −1 ).

Limitations
The applicability  throughout the ablation season, so that the effect of (i) and (ii) can be reduced. We use geodetic estimates for model calibration to limit the effect of (iv), (v) and (vi). Summer accumulation and all year round ablation are not taken in account for calibration, but they are included in the modelled daily mass balance. Therefore, it allows quantifying the annual mass balance for (iii) but makes the interpretation of seasonal to sub-seasonal mass balance time series more uncertain.

Debris cover
The influence of debris cover on glacier mass balance is complex and obscures the climate-related response of glacier mass changes. The approach applied here does not ac-

Glacier outlines
The glacier mass balance depends upon both, the climate and the configuration of the glacier geometry, each of which is a continuous function of time (Elsberg et al., 2001).
The so-called "reference" mass balance removes the surface change effects by holding the glacier surface constant through time, and thus is the climatically relevant mass balance.
However, it does not precisely quantify the actual annual mass gain and loss for a nonupdated outline. The so-called "conventional" mass balance takes in account the glacier geometry change and is thus relevant for quantifying meltwater release during a mass balance year (Harrison et al., 2005). Here, we chose constant glacier outlines taken from :

X -15
Asia date from the onset of the century and therefore give a good basis to calculate the reference mass balance from 1999/00 to 2017/18. Some regions (Eastern Pamir, Eastern Tien Shan, Dzungarsky Alatau and the Eastern part of the Central Tien Shan) contain outlines not only from the beginning of the study period but also from 2007 and 2009.
The choice of the reference surface affects the reference mass balance. It has been shown that the trend in the difference between the glacier-wide cumulative and the reference mass balance will grow approximately linearly with time but the difference between the glacier-wide annual balances on two different reference surfaces will be approximately constant over the years (Elsberg et al., 2001). For South Cascade glaciers, a difference in mass balance of 16% for the period 1970 to 1997 between conventional and reference balance was found (Elsberg et al., 2001). Following evolution, the difference in the mass balances using outlines from 2000 and 2009, respectively, would result in a difference of mass balance of less then 5%. The cumulative mass balance however is more sensitive to the choice of the reference area (Elsberg et al., 2001).

Model calibration
The model calibration is based on the information content of the transient snowline observation. In order to apply the approach on glaciers without direct measurements, the use of additional ground-based information was avoided. Thus, the two calibrated parameters are not strictly independent; meaning that an overestimated accumulation rate can be compensated by a too high ablation rate without degrading the agreement with the transient snowline observations. Erroneous seasonal components do not affect X -16 : excessively the calculated annual balance as they equal out but introduce considerable uncertainties for the seasonal to sub-seasonal model outputs. An attempt to overcome this problem was made by iterative calibration, where both parameters were adjusted in parallel (Barandun et al., 2018). With this iterative calibration, the model tends towards an optimal solution for both parameters without the need of additional observations. Because winter measurements are sparse for the region a detailed validation of the seasonal components was so far not possible, however, for the few years and few glaciers for which winter measurements are available, the calibration procedures provided satisfying results (Barandun et al., 2018). Nonetheless, at the current stage of the work, interpretation of seasonal to sub-seasonal mass balance time series has largely been avoided and more work is needed to properly evaluate the model performance on increased temporal resolution.    Figure S6.

SECOND-ORDER ADJUSTMENT FIRST-ORDER MODEL CALLIBRATION
Calibration procedure to obtain an ideal combination of DDF snow and C prec . In a first step, DDF snow and C prec are optimised through comparison to snowline observations until a good solution for both parameters was found. For the initial value of DDF snow , the best value of C prec was determined constraining the modelled cumulative melt M elt modelled at the snowline position to agree with the modelled winter snow accumulation Accu modelled for the same location. Then, the model performance was evaluated to minimise RM SE SCAF of the modelled and observed snow-covered area fractions (SCAF). This is repeated until no improvement in RMSE is obtained. In a second step, the modelled mass balance from 2004 to 2012 is compared to the multi-annual geodetic mass balance, and the relation between the degree day factor of snow and ice (R DDF ) and the precipitation gradient δP/δz were readjusted until the results of both approaches agree within the uncertainty range of the geodetic estimate.  Figure S12.
Year-to-year variability by means of standard deviation ( Figure S13. In order to investigate possible relationships of the mass balance with the key meteorological drivers, precipitation and temperature, we used the independent dataset High Asia Refined analysis (HAR) version 1.4 (Maussion et al., 2014). A. Monthly precipitation from HAR30 Reanalysis for the two catchments, for a below-average (red), and an above-average (red) mass balance year. B. Same but for monthly temperature.
X -34 :   Table S1. Comparison of modelled mass balances with geodetic estimates of other studies (Shean et al., 2020;Hugonnet et al., 2021).     Table S4. Monthly discharge (m 3 s −1 ) of the Gunt and Naryn Rivers as average, and for extremely negative and positive mass balance years obtained through Tajik Hydromet Office and Central-Asian Institute for Applied Geosciences. Baseline for the calculation of monthly means is the period 1999/00 to 2012/13 (Gunt) and 1999/00 to 2017/18 (Naryn), respectively. : X -37 Table S5. Total area, amount and median elevation of all glaciers, and number of