Spatial Heterogeneity in Glacier Mass-Balance Sensitivity across High Mountain Asia

Mass balance of glaciers in High Mountain Asia (HMA) varies substantially across the region. While the spatial variability is attributed to differences in climatic setting and sensitivity of these glaciers to climate change, an assessment of these factors to date has only been performed on a small sample of glaciers and a small set of climate perturbation scenarios. To advance the assessment to larger datasets, we first reconstruct the time series of reference-surface mass balance for 1952–2014 periods using an empirical model calibrated with observed mass balance from 45 glaciers across the HMA. Forcing the model with a set of independent stepwise changes of temperature (±0.5 K to ±6 K) and precipitation (±5% to ±30%), we assess the reference-surface mass balance sensitivity of each glacier in the sample. While the relationship between the change in mass balance and the change in precipitation is linear, the relationship with the change in temperature is non-linear. Spatial heterogeneity in the simulated mass balance sensitivities is attributed to differences in climatic setting, elevation, and the sensitivity of mass-balance profile (gradient) to changes in temperature and precipitation. While maritime and low-lying continental glaciers show high sensitivity to temperature changes and display a uniform mass-balance sensitivity with elevation, the high-lying continental glaciers show high sensitivity to precipitation changes and display a non-uniform mass-balance sensitivity with elevation. Our analysis reveals the dominant drivers of spatial variability in the mass balance sensitivity across the region: temperature as a single driver for maritime glaciers, and a superposition of temperature, precipitation seasonality, and snow/rain differentiation for continental glaciers. Finally, a set of sensitivity tests with perturbed model parameters confirms the robustness of our results. The model’s ability and robustness to resolve spatial patterns in the sensitivities and their drivers implies that simple modeling approaches remain a powerful tool for analyzing glacier response to climate change in HMA.


Introduction
Over the last few decades, rates of climate warming in High Mountain Asia (HMA) are among the highest on Earth [1].As high-elevation environments warm faster than low elevation, the retreat of mountain glaciers has been accelerating [2].Increased glacier melting during the twentieth century has

Study Area
HMA covers a large fraction of countries in Asia, spanning from Altai and Tianshan Mountains to the Tibetan Plateau (TBP) and surroundings (Figure 1), and includes parts of North Asia, the entirety of Central Asia, and Western and Eastern South Asia.The overall area is 9.86 × 10 4 km 2 , which is approximately 20% of the global area of mountain glaciers (excluding ice sheets of Greenland and Antarctica) [26].
Atmospheric circulation patterns over HMA are characterized by westerlies during winters and the Indian monsoon during summers.These synoptic patterns exert substantial effects on glaciers through the annual distribution of precipitation and temperature, therefore controlling the glacier type (summer-versus winter-accumulation) and spatial distribution of glaciers across the region [8].The Siberian anticyclonic circulation also influences glaciers in the Altai region and part of the Tianshan regions [27].Glaciers in the Qilian Mountains and the eastern TBP are influenced by the East Asia monsoon, while the continental climate dominates the interior of TBP.Interaction between the local climate and the South Asian monsoon drives the Karakoram anomaly [28].
According to the glacier temperature, glaciers in HMA can be categorized into three types: maritime glaciers, subcontinental glaciers and continental glaciers (Figure 1) [29].To systematically assess glacier mass balance sensitivity, 45 monitored glaciers are selected from different atmospheric circulation regions and different climatic setting (Figure 1).These 45 glaciers are unevenly distributed within the HMA: the majority of glaciers are located in the Tianshan Mountains and the Central Himalaya Mountains, while fewer are located in the inner and eastern TBP.Except for the West Kunlun and the interior of TBP, at least one of these glaciers is located in each climate setting.More detailed information on the 45 glaciers is provided in the Supplementary Table S1.
Water 2018, 10, x FOR PEER REVIEW 3 of 21 which is approximately 20% of the global area of mountain glaciers (excluding ice sheets of Greenland and Antarctica) [26].Atmospheric circulation patterns over HMA are characterized by westerlies during winters and the Indian monsoon during summers.These synoptic patterns exert substantial effects on glaciers through the annual distribution of precipitation and temperature, therefore controlling the glacier type (summer-versus winter-accumulation) and spatial distribution of glaciers across the region [8].The Siberian anticyclonic circulation also influences glaciers in the Altai region and part of the Tianshan regions [27].Glaciers in the Qilian Mountains and the eastern TBP are influenced by the East Asia monsoon, while the continental climate dominates the interior of TBP.Interaction between the local climate and the South Asian monsoon drives the Karakoram anomaly [28].
According to the glacier temperature, glaciers in HMA can be categorized into three types: maritime glaciers, subcontinental glaciers and continental glaciers (Figure 1) [29].To systematically assess glacier mass balance sensitivity, 45 monitored glaciers are selected from different atmospheric circulation regions and different climatic setting (Figure 1).These 45 glaciers are unevenly distributed within the HMA: the majority of glaciers are located in the Tianshan Mountains and the Central Himalaya Mountains, while fewer are located in the inner and eastern TBP.Except for the West Kunlun and the interior of TBP, at least one of these glaciers is located in each climate setting.More detailed information on the 45 glaciers is provided in the Supplementary Table 1.

Model Driving Data
To reconstruct the mass balance time series of the selected 45 glaciers, we used temperature and precipitation data from Climatic Research Unit (CRU).CRU TS v. 4.01 is a global-scale gridded climate dataset based on the station data interpolation with a spline-fitting technique [30].We used

Model Driving Data
To reconstruct the mass balance time series of the selected 45 glaciers, we used temperature and precipitation data from Climatic Research Unit (CRU).CRU TS v. 4.01 is a global-scale gridded climate dataset based on the station data interpolation with a spline-fitting technique [30].We used the monthly data provided for 1901-2016 at the spatial resolution of 0.5 • × 0.5 • .An earlier study has shown an agreement between CRU precipitation and rain gauge measurements from 29 meteorological stations in the high-elevation TBP [31].Despite underestimating the annual totals, CRU data has been shown to capture the inter-annual precipitation variability and spatio-temporal patterns within the region [31].Furthermore, previous studies successfully used CRU dataset to model glacier mass balance in HMA and other regions worldwide [32][33][34].CRU also provides one of the longest global climatology reconstructions for the 20th century; global reanalysis datasets, for example, ERA-40 and ERA-Interim do not reach further back than late 1950s) [35,36].Nevertheless, our initial modeling results showed poor performance with the use of CRU data at four of our 45 glaciers: Urumqi No. 1, Muztag Ata, Naimona'nyi, and Baishui Glacier No. 1.For these four glaciers, we replaced the CRU data with temperature and precipitation from the national weather stations in the glaciers' vicinity (Figure 1): Daxigou station (3539 meters above sea level (m a.s.l.), 86.5 Surface area for all glaciers within China is provided by the second Chinese Glacier Inventory (CGI) which is an updated version of the first CGI [37,38].For most glaciers in CGI, the area is derived from the 2006-2010 period.For the glaciers outside China we used the Randolph Glacier Inventory (RGI) 5.0, which is based on the satellite imagery from 1999-2010 [26].Glacier hypsometry (area versus elevation) was calculated for 50 m elevation intervals using the Shuttle Radar Topography Mission (SRTM) DEM [39].
The World Glacier Monitoring Service (WGMS) collects and publishes observed mass balance data for glaciers worldwide [40].At present, WGMS provides mass balance data for 65 glaciers in HMA.For our study, we selected the glaciers with at least three years of mass balance record.After initial model runs, we excluded a set of glaciers with relatively large errors between the modeled and observed mass balance time-series.The errors could not be explained by poorly constrained model parameters or poor performance in CRU data.The eliminations led to the selection of 45 glaciers in total (Figure 1), whose length of observed record varies from three to 54 years.For example, Urumqi No. 1 glacier has the longest record (27 years) among the Chinese glaciers, followed by glaciers in the Altai and Tianshan (average record of 18 years).More information on the mass balance record is provided in the Supplementary Table S1.Since WGMS has an incomplete mass balance record for Siachen glacier, we used the geodetic mass balance observations derived from Cartosat-1 and SRTM DEM [41].

Mass Balance Model
To reconstruct mass balance time series for each of our 45 selected glaciers, we adopted the model from Radic and Hock [25,42].Key input data are monthly precipitation and temperature.Area-weighted specific mass balance (B) for the whole glacier is computed as a sum of specific mass balance (b) of each elevation band (i): where S i is surface area of an elevation band, and n is the total number of elevation bands with elevation interval of 50 m.For each elevation band, we calculated monthly specific mass balance (mm water equivalent) as: where a i is surface ablation (negative), c i is mass accumulation (positive), and R i refers to melt water refreezing (positive).The ablation (mm w.e.) of snow and ice is calculated based on a positive degree-day model [42]: where f snow/ice is a degree-day factor for snow and ice (mm w.e.d −1• C −1 ), respectively, T i is the near-surface air temperature ( • C), as a monthly mean, and dt converts to a one-day time step in a discretized representation of the integral.
Monthly accumulation for each elevation band, c i (mm w.e.) is calculated as: where T snow represents the threshold temperature to differential between snow and rainfall, and P i is monthly precipitation.If T i is below T snow , a solid precipitation is assumed contributing to mass accumulation, otherwise rainfall is assumed and the accumulation is set to zero (δ m is constant).The annual potential refreezing depth for each elevation band, R i (cm w.e.), is calculated though the empirical relationship with annual mean air temperature, T a ( • C), as [43]: where the lower boundary of annual potential depth over the whole glacier is 0, while an upper boundary is the snow melt water (cm w.e.).Monthly melt water does not flow away until accumulated melt in a mass balance year exceeds the annual potential depth.
To interpolate CRU monthly temperature to each elevation band, we utilized two temperature lapse rates, Tlap and Glap.The Tlap corresponds to the 'statistical lapse rate', i.e., the temperature correction between the elevation of the CRU grid cell (h CRU ) and the maximum elevation band (h max ) of the glacier.Glap represents the temperature lapse rate along the glacier surface.The monthly air temperature T i at each elevation band of height h is calculated as: To upscale the CRU monthly precipitation from the grid cell to the maximum elevation band (h max ), we used a precipitation correction factor, kp, while a precipitation gradient, dpre, is used to interpolate precipitation along the glacier surface.The dpre is assigned as a percentage of precipitation decrease with every 100 m drop and is calculated from the top to the bottom of each glacier.Thus, the precipitation for each elevation band is computed as: The glacier mass balance sensitivities to temperature change (∆B/∆T) and precipitation change (∆B/∆P) are calculated as: ∆B Water 2019, 11, 776 where B is the modeled mean annual mass balance (m w.e.year −1 ) for the period of 1952-2014, B(±∆T) is the same B but in response to a step-wise temperature change in the range of −6 K to +6 K with increment of 0.5 K, and B(±∆P) is same B, but in response to a step-wise precipitation change of −30% P to +30% P with an increment of 5% P. ∆T and ∆P are introduced as a uniform stepwise change to the original monthly temperature and precipitation time series throughout the 63-year period.Annual mass balance is derived for a hydrological year, i.e., from October to September of the following year.Summer mass balance is calculated from April to September, and winter mass balance is from October to March of the following year.

Model Parameter Optimization
We follow the model optimization procedure from Radic and Hock [25] to obtain the best match between modeled and observed summer (winter) mass balance time series.Note that, in the original study [25], the glacier area is kept fixed during the calibration period, while the volume-area-length scaling method, accounting for area and volume evolution, is applied only in the projection period.For the consistency with the original study, we kept the area of each glacier fixed throughout the 63-year reconstruction period.Therefore, our reconstructed mass balance time series represent the reference-surface mass balance [44].We set the initial range for each parameter following the original study [25] and then applied an optimization method that minimizes the root-mean-square-error (RMSE) between the modeled results and the observed seasonal (winter and summer) mass balance time series for each glacier separately.A summary of the statistics (min, max, mean, and standard deviation) across the 45 glaciers for each parameter is shown in Table 1, while more details are provided in the Supplementary Table S2.Among the 45 glaciers, the largest RMSE between modeled and observed annual mass balance time series is for Parlung No. 390 glacier (0.63 m w.e year −1 ) located in the Hengduan Mountains.In general, the shorter the observed mass balance record means the larger the RMSE.The smallest RMSE (0.06 m w.e.year −1 ) is found for Zhadang glacier in the Nyainqêntanglha.The average RMSE across 45 glaciers is 0.31 m w.e year −1 for annual mass balance (Supplementary Table S2), 0.42 m w.e. for summer, and 0.38 m w.e. for winter mass balance (Figure 2).A more detailed comparison between modeled and observed annual mass balance time series for each glacier is shown in Supplementary Figure S2.The poorest model performance is observed for Neh Nar, Yala, Pokalde, and Changmekhangpu glaciers in the Himalaya.We attribute this poor performance to the inability of CRU data to capture the complex precipitation patterns in the mountainous region, and to the uncertainty in extrapolating from the in-situ mass balance measurements to the whole glacier.We note that all available years of mass balance observations are used in the model optimization.A more statistically robust way to perform the optimization is to use, for example, cross-validation methods applied over an independent sample of observations from the optimization sample.These methods, however, rely on having a larger dataset than available here.

Mass Balance Sensitivity
The mass balance sensitivities of each glacier to temperature changes of −6 K to 6 K and precipitation changes of −30% to 30% are shown in Figure 3.The results show a difference in the patterns between sensitivity to temperature and sensitivity to precipitation: while the change in mass balance linearly increases with the increase in precipitation perturbation, the relationship between the change in mass balance and the change in temperature perturbation is non-linearly.According to the temperature sensitivity, maritime glaciers are the most sensitive and continental glaciers are the least sensitive.Note that once the temperature perturbation drops below −3 K, the change in mass balance stays unaltered for most of the glaciers (Figure 3a).In other words, a uniform temperature drop of more than −3 K leads to below-freezing conditions year-round: all annual precipitation falls as snow and there is no ablation.Thus, any further decrease in temperature would not alter glacier mass balance.On the other hand, for a temperature perturbation of +0.5 K, the absolute change of mass balance sensitivity increases linearly with the temperature perturbation (Figure 3b).We find that the average temperature sensitivity across those 45 glaciers is −0.89 ± 0.36, −1.21 ± 0.46, and −1.53 ± 0.56 m w.e.year −1 under temperature increase +1.5 K, +2 K, and +2.5 K, respectively.The magnitude of precipitation sensitivity shows a linear dependency on precipitation change (Figure 3c), which translates into a constant sensitivity per % of precipitation change per year.The average precipitation sensitivity across the 45 glaciers is −0.23 and 0.22 m w.e.year −1 , for −30% and 30% precipitation change, respectively (Figure 3d).

Mass Balance Sensitivity
The mass balance sensitivities of each glacier to temperature changes of −6 K to 6 K and precipitation changes of −30% to 30% are shown in Figure 3.The results show a difference in the patterns between sensitivity to temperature and sensitivity to precipitation: while the change in mass balance linearly increases with the increase in precipitation perturbation, the relationship between the change in mass balance and the change in temperature perturbation is non-linearly.According to the temperature sensitivity, maritime glaciers are the most sensitive and continental glaciers are the least sensitive.Note that once the temperature perturbation drops below −3 K, the change in mass balance stays unaltered for most of the glaciers (Figure 3a).In other words, a uniform temperature drop of more than −3 K leads to below-freezing conditions year-round: all annual precipitation falls as snow and there is no ablation.Thus, any further decrease in temperature would not alter glacier mass balance.On the other hand, for a temperature perturbation of +0.5 K, the absolute change of mass balance sensitivity increases linearly with the temperature perturbation (Figure 3b).We find that the average temperature sensitivity across those 45 glaciers is −0.89 ± 0.36, −1.21 ± 0.46, and −1.53 ± 0.56 m w.e.year −1 under temperature increase +1.5 K, +2 K, and +2.5 K, respectively.The magnitude of precipitation sensitivity shows a linear dependency on precipitation change (Figure 3c), which translates into a constant sensitivity per % of precipitation change per year.The average precipitation sensitivity across the 45 glaciers is −0.23 and 0.22 m w.e.year −1 , for −30% and 30% precipitation change, respectively (Figure 3d).To analyze the spatial heterogeneity in the mass balance sensitivity across the region, we use the scenarios with ±1 K temperature change and ±10% P precipitation change (Figure 4).Temperature sensitivities to −1 K range from 0.17 to 1.54 m w.e.year −1 K −1 , while for +1 K they range from −0.23 to −1.43 m w.e.year −1 K −1 .For temperature sensitivities to +1 K, our average value of −0.42 m w.e.year −1 K −1 is consistent with the average of −0.56 and −0.49m w.e.year −1 K −1 found for Altai region [25,45].For Tianshan, our average value of −0.47 m w.e.year −1 K −1 is slightly higher than −0.42 m w.e.year −1 K −1 previously found [25], while for Himalaya, our average value of −0.62 m w.e.year −1 K −1 is slightly smaller than −0.68 m w.e.year −1 K −1 [25].On individual glacier scale, our results agree with previously assessed sensitivities for Shumskiy and Urumqi No. 1 glacier, but are smaller compared to the assessments for other glaciers [17].The average temperature sensitivity of −0.59 m w.e.year −1 K −1 for the 45 glaciers is lower than the mean sensitivity found for the glaciers in Scandinavia (−0.61 m w.e.year −1 K −1 ) [46] and New Zealand (−1.1 m w.e.year −1 K −1 ) [11], but higherthan the global mean sensitivity (−0.41 m w.e.year −1 K −1 ) [47].Spatial heterogeneity of the mass balance sensitivity (Figure 4) corroborates previous findings that maritime glaciers are most sensitive to temperature and precipitation changes, continental glaciers are least sensitive, and subcontinental glaciers display a more complex pattern [16].To analyze the spatial heterogeneity in the mass balance sensitivity across the region, we use the scenarios with ±1 K temperature change and ±10% P precipitation change (Figure 4).Temperature sensitivities to −1 K range from 0.17 to 1.54 m w.e.year −1 K −1 , while for +1 K they range from −0.23 to −1.43 m w.e.year −1 K −1 .For temperature sensitivities to +1 K, our average value of −0.42 m w.e.year −1 K −1 is consistent with the average of −0.56 and −0.49m w.e.year −1 K −1 found for Altai region [25,45].For Tianshan, our average value of −0.47 m w.e.year −1 K −1 is slightly higher than −0.42 m w.e.year −1 K −1 previously found [25], while for Himalaya, our average value of −0.62 m w.e.year −1 K −1 is slightly smaller than −0.68 m w.e.year −1 K −1 [25].On individual glacier scale, our results agree with previously assessed sensitivities for Shumskiy and Urumqi No. 1 glacier, but are smaller compared to the assessments for other glaciers [17].The average temperature sensitivity of −0.59 m w.e.year −1 K −1 for the 45 glaciers is lower than the mean sensitivity found for the glaciers in Scandinavia (−0.61 m w.e.year −1 K −1 ) [46] and New Zealand (−1.1 m w.e.year −1 K −1 ) [11], but higherthan the global mean sensitivity (−0.41 m w.e.year −1 K −1 ) [47].Spatial heterogeneity of the mass balance sensitivity (Figure 4) corroborates previous findings that maritime glaciers are most sensitive to temperature and precipitation changes, continental glaciers are least sensitive, and subcontinental glaciers display a more complex pattern [16].

Sensitivity of Mass Balance Gradients
To further analyze the differences among the mass balance sensitivities of individual glaciers, we selected nine glaciers of different glacier type and from different subregions of HMA (Figure 5).We find that temperature sensitivity for the three maritime glaciers (Yala, Parlung No. 10, and Baishui Glacier No. 1) is roughly uniform with elevation, i.e., mass balance gradient sensitivity is linear.On the other hand, the sensitivity for the two continental glaciers (Muztag Ata and Siachen) shows non-uniform dependence with elevation, corroborating previous findings [24].Our findings reveal that the maritime glaciers lose their accumulation areas as the temperature increases, thus transforming the original accumulation-dominated elevation band into ablation-dominated bands.This process, as the temperature further increases, will turn the entire glacier into the ablation area, shifting the ELA above the top of the glacier (Figure 6).By contract, the accumulation-dominated elevation bands of the continental glaciers are less sensitive, leading to a smaller shift in the ELA (Figure 6).This gradient sensitivity and ELA behavior demonstrate that the summer-accumulation type maritime glaciers are controlled by temperature, whereas the continental glaciers are controlled by temperature and precipitation.Mass balance at low-elevation is mainly controlled by temperature, while at high-elevations by precipitation.For the precipitation sensitivity, the mass balance at high-elevations is more sensitive to precipitation than at low elevations, once again suggesting that the high-elevation bands are mainly affected by precipitation.

Sensitivity of Mass Balance Gradients
To further analyze the differences among the mass balance sensitivities of individual glaciers, we selected nine glaciers of different glacier type and from different subregions of HMA (Figure 5).We find that temperature sensitivity for the three maritime glaciers (Yala, Parlung No. 10, and Baishui Glacier No. 1) is roughly uniform with elevation, i.e., mass balance gradient sensitivity is linear.On the other hand, the sensitivity for the two continental glaciers (Muztag Ata and Siachen) shows non-uniform dependence with elevation, corroborating previous findings [24].Our findings reveal that the maritime glaciers lose their accumulation areas as the temperature increases, thus transforming the original accumulation-dominated elevation band into ablation-dominated bands.This process, as the temperature further increases, will turn the entire glacier into the ablation area, shifting the ELA above the top of the glacier (Figure 6).By contract, the accumulation-dominated elevation bands of the continental glaciers are less sensitive, leading to a smaller shift in the ELA (Figure 6).This gradient sensitivity and ELA behavior demonstrate that the summer-accumulation type maritime glaciers are controlled by temperature, whereas the continental glaciers are controlled by temperature and precipitation.Mass balance at low-elevation is mainly controlled by temperature, while at high-elevations by precipitation.For the precipitation sensitivity, the mass balance at high-elevations is more sensitive to precipitation than at low elevations, once again suggesting that the high-elevation bands are mainly affected by precipitation.

Model Sensitivity
Main sources of uncertainty in the estimates of mass balance sensitivity are rooted in the choice of model optimization and the choice of parameter values.While the former depends on the accuracy of input (temperature and precipitation) and output (observed mass balance) data in the model [48], the latter can be tested by perturbing each parameter value and assessing the model sensitivity to these perturbations.
Here, we run the model for each glacier with perturbed value of each parameter separately.We chose to perturb each parameter by ±20% of the original (optimized) value, while the values of all other parameters are kept fixed at their optimized values.The model sensitivity to parameter perturbations is calculated as the difference (RMSE and bias) between the modeled and observed mass balance time series averaged across the 45 glaciers (Table 2).We find that the model is most sensitive to the temperature lapse rate (Tlap), with perturbation of −20% resulting in the model sensitivity of −1.59 m w.e.year −1 , while the model is least sensitive to snow/rain threshold temperature (Tsnow), with the perturbation of −20% yielding −0.45 m w.e.year −1 (Table 2).The results also reveal that different parameter combinations, within 20% difference from the optimized values, yield a relatively good match (small RMSE) to the observed mass balance.As expected, the smallest average RMSE across the 45 glaciers is obtained with the optimized (original) parameter values.
Next we assessed the changes in the spatial variability of mass balance sensitivity when the model is run with the perturbed parameter values.The model is run separately with the following perturbations in each parameter: −20% (dpre), 20% (Glap), −20% (Tlap), −20% (Tsnow), −20% (fsnow), −20% (fice), and −20% (kp).Note that we do not consider the runs with perturbed Tlap due to substantially higher model sensitivity to this parameter relative to the other parameters.The mass balance sensitivity is assessed for the temperature change of +1 K and the precipitation change of +10% for each of the 45 glaciers (Figure 7).Results show that most parameters have a strong influence on the mass balance sensitivity, with the strongest influence by temperature lapse rate (Glap), yielding up to 0.9 m w.e.year −1 difference the new sensitivity and original sensitivity value (Table 2).While the parameter perturbations led to some large differences in mass balance

Model Sensitivity
Main sources of uncertainty in the estimates of mass balance sensitivity are rooted in the choice of model optimization and the choice of parameter values.While the former depends on the accuracy of input (temperature and precipitation) and output (observed mass balance) data in the model [48], the latter can be tested by perturbing each parameter value and assessing the model sensitivity to these perturbations.
Here, we run the model for each glacier with perturbed value of each parameter separately.We chose to perturb each parameter by ±20% of the original (optimized) value, while the values of all other parameters are kept fixed at their optimized values.The model sensitivity to parameter perturbations is calculated as the difference (RMSE and bias) between the modeled and observed mass balance time series averaged across the 45 glaciers (Table 2).We find that the model is most sensitive to the temperature lapse rate (Tlap), with perturbation of −20% resulting in the model sensitivity of −1.59 m w.e.year −1 , while the model is least sensitive to snow/rain threshold temperature (Tsnow), with the perturbation of −20% yielding −0.45 m w.e.year −1 (Table 2).The results also reveal that different parameter combinations, within 20% difference from the optimized values, yield a relatively good match (small RMSE) to the observed mass balance.As expected, the smallest average RMSE across the 45 glaciers is obtained with the optimized (original) parameter values.
Next we assessed the changes in the spatial variability of mass balance sensitivity when the model is run with the perturbed parameter values.The model is run separately with the following perturbations in each parameter: −20% (dpre), 20% (Glap), −20% (Tlap), −20% (Tsnow), −20% (fsnow), −20% (fice), and −20% (kp).Note that we do not consider the runs with perturbed Tlap due to substantially higher model sensitivity to this parameter relative to the other parameters.The mass balance sensitivity is assessed for the temperature change of +1 K and the precipitation change of +10% for each of the 45 glaciers (Figure 7).Results show that most parameters have a strong influence on the mass balance sensitivity, with the strongest influence by temperature lapse rate (Glap), yielding up to 0.9 m w.e.year −1 difference the new sensitivity and original sensitivity value (Table 2).While the parameter perturbations led to some large differences in mass balance sensitivity on the scale sensitivity on the scale of individual glaciers, the overall pattern of sensitivity, with a roughly identical trend of change, is preserved in all the runs (Figure 7).

Drivers of Mass Balance Sensitivity
Previous studies indicated that mass balance sensitivity is closely related to mass balance amplitude [11,49], precipitation [12], and continentality index (CI, i.e., the temperature difference Water 2019, 11, 776 13 of 21 between the highest and lowest monthly temperature within a year) [13,50].Precipitation has an important influence on mass balance sensitivity through the mass turnover [12].To some extent, precipitation affects the mass balance sensitivity through changes in the ratio of liquid versus solid precipitation, which is again controlled by temperature.As temperature decreases, summer snowfall on a glacier surface will increase and change the mass balance.Thus, the mass balance sensitivity to temperature is significantly correlated with precipitation (Figure 8a).However, when the decrease in temperature exceeds −3 K, all precipitation becomes snowfall and any further temperature decrease does not alter the snowfall proportion in the total precipitation.Hence, the correlation coefficient cease to vary for ∆T < −3 K, while for ∆T > −3 K the correlation decreases because the effect of precipitation is reduced and glacier ablation dominates temperature sensitivity (Figure 8a,b).

Factors Influencing Temperature Sensitivity
Our results confirm that the degree-day factor and precipitation have significant effect on the mass balance sensitivity of glaciers in HMA.Furthermore, we showed that the ELA can influence the mass balance sensitivity by altering the accumulation rate and, particularly, the gradient sensitivity (Figure 5).The maritime glaciers have a relatively small accumulation area and a high ELA (high relative to the bottom of the glacier): as their near-surface temperature increases, the accumulation area shrinks and the ELA can shift above the top of the glacier.The linearity in the mass-balance gradient with elevation, therefore, reflects the sensitivity of the ablation area to temperature change.By contrast, the continental and subcontinental glaciers have larger accumulation area and low ELA (low relative to the bottom of the glacier).The non-linearity in their mass-balance gradient with elevation, therefore, reflects the differences between sensitivities of ablation and accumulation area to temperature change.Hence, we conclude that the position of ELA relative to the top of the glacier has an important influence on the shape of the temperature sensitivity gradient.
As our results show, the degree-day factor is a key parameter influencing the variability of In general, the degree-day factor indicates the melting extent of glaciers and affects temperature sensitivity through the ablation process.Our results for ∆T > 0 K show that the ablation is a main driver of mass balance sensitivity to temperature (Figure 8c,d).Furthermore, for 0 < ∆T < 3 K, the correlation between the temperature sensitivity and the degree day factor slightly decreases due to the influence of snowfall.For ∆T > 3 K, the temperature sensitivity is mainly affected by the degree day factor.We conclude that precipitation and the degree day factor affect the mass balance sensitivity under temperature changes within ±3 • C.
It has been shown that CI influences temperature sensitivity because higher the summer temperatures, larger the glacier ablation.Furthermore, a temperature increase of +1 K that occurs around the melting temperature can significantly alter the ratio of liquid versus solid precipitation [11].
Water 2019, 11, 776 14 of 21 Applying the same correlation analysis to the mass-balance sensitivity and CI, reveals that CI is significantly correlated with temperature sensitivity for −1 < ∆T < 2 K (Figure 8e,f).Therefore, the maritime glaciers in the monsoon zones of HMA, with relatively low CI, have larger summer ablation than other glaciers in the region.
Terrain and glacier topography can also influences the mass balance sensitivity through factors such aspect, slope, and Aspect can drive the sensitivity by dictating the duration and amount of glacier surface exposure to solar radiation [51].The steepness of the slope and elevation range can dictate the temperature and precipitation gradients along the glacier surface [11].In our sample of 45 glaciers, however, we did not find any significant correlation between the balance sensitivity and these factors.

Factors Influencing Temperature Sensitivity
Our results confirm that the degree-day factor and precipitation have significant effect on the mass balance sensitivity of glaciers in HMA.Furthermore, we showed that the ELA can influence the mass balance sensitivity by altering the accumulation rate and, particularly, the gradient sensitivity (Figure 5).The maritime glaciers have a relatively small accumulation area and a high ELA (high relative to the bottom of the glacier): as their near-surface temperature increases, the accumulation area shrinks and the ELA can shift above the top of the glacier.The linearity in the mass-balance gradient with elevation, therefore, reflects the sensitivity of the ablation area to temperature change.By contrast, the continental and subcontinental glaciers have larger accumulation area and low ELA (low relative to the bottom of the glacier).The non-linearity in their mass-balance gradient with elevation, therefore, reflects the differences between sensitivities of ablation and accumulation area to temperature change.Hence, we conclude that the position of ELA relative to the top of the glacier has an important influence on the shape of the temperature sensitivity gradient.
As our results show, the degree-day factor is a key parameter influencing the variability of mass balance sensitivity in HMA.The glaciers with large degree-day factors are regarded as having a strong melt potential.Across our glacier sample, the maritime glaciers display the highest degree-day factors in the region, while the lowest factors are found for the continental glaciers.The degree-day factors for the subcontinental glaciers substantially vary within the region.The spatial pattern of the degree-day factors is consistent with the patterns found in previous studies [52], highlighting the heterogeneity in the glacier ablation potential across the region.In the maritime climatice setting under the influence of Indian monsoon, the degree-day factors are large and the glacier mass loss is dominated by melting, while in the arid and cold environment dominated by westerlies, the degree-day factors are small and the glacier mass loss is controlled by sublimation [18,53,54].In order to provide a physics-based explanation for the glacier-to-glacier variability in the degree-days factors, one would need to look into the surface energy balance of each glacier [53], which is beyond the scope of this study.
In general, sensitivity of maritime glaciers to climate change is shown to be linked to high annual precipitation and subsequent mass turnover [12,55].More important than the link with the annual precipitation, however, is the link with precipitation seasonality [15,56,57].Figure 9b shows the area coverage of the winter-accumulation type glaciers (June-September), for which summer precipitation accounts for less than 60% of the annual precipitation.These glaciers are located in the Tianshan, Pamir, Karakoram, and Western Kunlun.The snowfall proportion for the summer-accumulation type glaciers (for which summer precipitation accounts for more than 60% of the annual precipitation) is highly dependent on the temperature changes because their near-surface air temperatures are close to the rain/snow threshold temperature.Thus, a small temperature change can trigger a large change in the annual snowfall proportion at these glaciers [11,21,51].Under a temperature increase of +1 K, the largest reduction in snowfall was detected for the summer-accumulation type glaciers (Figure 9b): the reduction of roughly 50% in snowfall has reduced the glacier mass accumulation and decreased the surface albedo further intensifying glacier melting [12,23].Furthermore, a temperature increase of +1 K enhanced the duration of annual ablation by one month on average.On the other hand, no obvious reduction in snowfall proportion to total precipitation nor a reduction in the ablation duration was detected for the winter-accumulation type glaciers under the same +1 K scenario.

Regional Patterns of Mass Balance And Sensitivity
Our sample of 45 glaciers with in situ mass balance observations has a bias toward small size glaciers.Thus, one needs to be careful when extrapolating the findings to all glaciers in the region, especially for the central Himalaya, where many glaciers with larger elevation band distributions reside [60].Overall, the pattern of spatial heterogeneity in glacier mass balance detected in our glacier sample agrees with the results based on ICEsat [5], GRACE [7], and ASTER [6].In particular, the patterns of mass gain for glaciers in the Karakoram and mass loss for glaciers in Himalaya and Hengduanshan are detected in all studies.In Karakoram, the glaciers mass change results in slightly accelerated glacier flow due to changes in gravitational driving stress controlled by changes in ice thickness, whereas glaciers in Himalaya and Hengduanshan show sustained slowdown due to ice thickness thinning, which can modify the glacier mass balance sensitivity to some extent [61].In addition, the winter snowfall increase and summer temperature decrease favor mass gain in Karakoram [62], while summer temperature increase and precipitation decrease drive the substantial mass loss in Himalaya and Hengduanshan [63].Our findings from the sample of 45 glaciers reveal that the key reasons for the positive mass balance in Karakoram and West Kunlun is the lower mass-balance sensitivity to temperature changes as their ELA sits relatively low.On the other hand, the mass balance of glaciers in Himalayas is very sensitive to temperature change, as their ELA sits relatively high.Figure 9c shows the spatially distributed differences in the refreezing proportion in response to +1 K scenario.Most refreezing occurs in the arid and cold environments, occupied by the continental glaciers, such as Muztag Ata and Siachen.The refreezing of glacier melt water in the firn layer inhibits the glacier mass loss.In the humid and warm environments, which are home to the maritime glaciers, the refreezing decreases as the temperature increases, which promotes further glacial melting.As a result, the maritime glaciers can become more sensitive as the temperatures keep rising.
As shown in Equation ( 2), the modeled mass balance is a superposition of three processes: surface ablation, accumulation and refreezing.We examined the change in each process under a scenario of +1 K warming (Figure 9d) and its impact on the mass balance sensitivity.Glacier ablation is clearly the main factor influencing temperature sensitivity and is controlled by the degree-day factor.The contribution of accumulation to temperature sensitivity, through the snowfall reduction, varies greatly across the region (from 0.8-46.2%)with the largest influence on the summer-accumulation type glaciers.The effect is particularly important for maritime glaciers, which are highly dependent on summer precipitation and receive the largest amounts of precipitation in the region (1453-3000 mm year −1 ) [58,59].The largest contribution of rising-temperature-induced snowfall reduction to temperature sensitivity is in Nyainqêntanglha and Himalaya (approximately 46.2%).In contrast, the majority of continental glaciers, exposed to the winter accumulation under the westerlies, have relatively low temperature sensitivity due to the rising-temperature-induced snowfall reduction.In conclusion, the spatial variability in the temperature sensitivity across the HMA can be largely explained by the variability in the degree-day factors ('glacier response') and the precipitation seasonality ('climate setting').

Regional Patterns of Mass Balance and Sensitivity
Our sample of 45 glaciers with in situ mass balance observations a bias toward small size glaciers.Thus, one needs to be careful when extrapolating the findings to all glaciers in the region, especially for the central Himalaya, where many glaciers with larger elevation band distributions reside [60].Overall, the pattern of spatial heterogeneity in glacier mass balance detected in our glacier sample agrees with the results based on ICEsat [5], GRACE [7], and ASTER [6].In particular, the patterns of mass gain for glaciers in the Karakoram and mass loss for glaciers in Himalaya and Hengduanshan are detected in all studies.In Karakoram, the glaciers mass change results in slightly accelerated glacier flow due to changes in gravitational driving stress controlled by changes in ice thickness, whereas glaciers in Himalaya and Hengduanshan show sustained slowdown due to ice thickness thinning, which can modify the glacier mass balance sensitivity to some extent [61].In addition, the winter snowfall increase and summer temperature decrease favor mass gain in Karakoram [62], while summer temperature increase and precipitation decrease drive the substantial mass loss in Himalaya and Hengduanshan [63].Our findings from the sample of 45 glaciers reveal that the key reasons for the positive mass balance in Karakoram and West Kunlun is the lower mass-balance sensitivity to temperature changes as their ELA sits relatively low.On the other hand, the mass balance of glaciers in Himalayas is very sensitive to temperature change, as their ELA sits relatively high.
We note that our mass balance model is calibrated locally, i.e., for each individual glacier separately, using seasonal (winter and summer) mass balance observations.Thus, each model parameter is glacier specific, which differs from the previous approach [42] where one of the parameters (temperature lapse rate) was tuned regionally to match modeled to observed regional mass balance time series.We showed that our locally-calibrated modeling approach reveals regional patterns in glacier mass balance and mass balance sensitivity that corroborate the assessments based on the remote sensing data [5][6][7].These findings, in addition to our findings on the non-linearity in glacier mass balance sensitivity across the region, advance our understanding of HMA glaciers.

Uncertainties in the Modeled Mass Balance and Sensitivity
For our mass balance modeling and sensitivity assessment, we kept the glacier area constant in time, thus operating with the 'reference-surface' mass balance.According to the glacier inventories used (CGI and RGI), the area uncertainty for each glacier is assumed to be 3.2-10% [37,64].Nevertheless, according to the CGI, glaciers in China have lost about 13% of their area during the period 1960s-2010 [65].The most pronounced glacier retreat and area reduction occurred in the Himalaya and Hengduan Mountains [8].The constant area and altitude distribution, as were assumed in our study, can yield an overestimated mass loss [66].Since the assumption neglects the area shrinkage, which mainly occurs at the lower-elevation bands that have the highest ablation rates along the glacier, our mass balance sensitivities are likely overestimated.
Despite the model's simplicity and the use of monthly time series of temperature and precipitation, the reconstructed mass balance time series of Muztag Ata glacier in Pamir and Zhangdang in the Nyainqêntanglha agree with the mass balance time series simulated with an energy balance [53,67] (see Supplementary Figures S1 and S2 for more details).In general, a temperature index model performs better on glaciers whose ablation is dominated by melt rather than sublimation.While ablation of glaciers in Pamir and Karakorum is dominated by sublimation, the agreement between our and the energy balance model gives additional confidence in our results.Nevertheless, our model does not take into account snow redistribution, topographic shading and variable albedo, which are the factors shown to be important for the mass balance modeling in this region [68].
Observations of annual mass balance, obtained by the in-situ glaciological method, are commonly regarded as the best data for model calibration [69].A typical measurement uncertainty in the annual mass balance for one site is considered to be 0.2 m w.e [70].Despite the increasing number of annual mass balance observations from geodetic and remote sensing methods for this region, we relied the model calibration with both summer and winter mass balance time series, which are more readily available from the in-situ measurements.Thus, we have limited our glacier sample in order to better constrain model parameters.A natural extension of this study would be to evaluate our results on glaciers with both in situ and remotely-observed mass balances.

Conclusions
The main goal of this study was to provide a systematic analysis of glacier mass balance sensitivity to temperature and precipitation changes across the HMA.To reach this goal, we first reconstructed the reference-surface mass balance time series for 1952-2014 period for 45 glaciers in the region.We used a simple mass balance model and, for each glacier, optimized the model parameters using the available in-situ observations of seasonal mass balances.To obtain the mass balance sensitivity of each glacier, we forced the model with a set of independent stepwise changes of temperature (±0.5 K to ±6 K) and precipitation (±5% to ±30%).The estimated mass balance sensitivity is shown to depend on the choice of parameter values in the model, especially to the choice of temperature lapse rate and precipitation correction factor.However, the key finding below are obtained in all model runs, i.e., with optimized and with perturbed parameter values: Mass balance sensitivity to temperature displays a substantially different pattern than the sensitivity to precipitation: while the change in reference-surface mass balance linearly increases with the change in precipitation, the relationship between the change in mass balance and the change in temperature is non-linear.
The temperature sensitivity displays substantial spatial heterogeneity across the region.The maritime glaciers are most sensitive to temperature change due to a rise in ELA above the uppermost glacier elevation, the continental glaciers are the least sensitive, and the subcontinental glaciers are in the middle.
The dependence of mass balance sensitivity on the climate regime reveals the dominant climatic drivers of the sensitivity: temperature as a single driver for the maritime glaciers, and a superposition of temperature, precipitation seasonality, and snow/rain differentiation for the continental glaciers.
Our study demonstrated the effectiveness of a simple empirical model of glacier mass balance, initially developed for the assessment of glacier contribution to global sea level rise, to quantify spatial heterogeneity in mass balance sensitivity across HMA.The model, originally designed for the global scale assessment, is shown to be successful in its application on the regional/local scale.Furthermore, despite the bias towards small glaciers in our glacier sample, the detected patterns of mass balance and its sensitivity are representative of those found for all glaciers in the region.These findings and the model's robustness demonstrated with the sensitivity tests imply that simple modeling approaches applied on glaciers with in-situ measurements still remain a powerful tool in studying complex patterns of glacier response to climate change across spatial and temporal scales.

Figure 1 .
Figure 1.Location of the 45 monitored glaciers in HMA used in this study (triangles) and four national weather stations (green stars).Polygons are enclosing the glaciers of a certain type: maritime type is within orange polygon, continental type within purple, while all outside are subcontinental types.Next to Mametova glacier are seven unmarked glaciers: Ts.Tuyuksu, Igly.Tuyuksu, Molodezhniy, Kara-Batkak, Mayakovskiy, Ordzhonikidze, and Visyachiy; next to Dunagiri glacier is one unmarked glacier: Naimona'nyi; and next to Pokakde glacier is one unmarked glacier: Changmekhangpu.

Figure 1 .
Figure 1.Location of the 45 monitored glaciers in HMA used in this study (triangles) and four national weather stations (green stars).Polygons are enclosing the glaciers of a certain type: maritime type is within orange polygon, continental type within purple, while all outside are subcontinental types.Next to Mametova glacier are seven unmarked glaciers: Ts.Tuyuksu, Igly.Tuyuksu, Molodezhniy, Kara-Batkak, Mayakovskiy, Ordzhonikidze, and Visyachiy; next to Dunagiri glacier is one unmarked glacier: Naimona'nyi; and next to Pokakde glacier is one unmarked glacier: Changmekhangpu.

Water 2018 ,
10, x FOR PEER REVIEW 7 of 21

Figure 2 .
Figure 2. Comparison between modeled and observed annual, summer, and winter balance for each glacier and each year in our dataset.

Figure 2 .
Figure 2. Comparison between modeled and observed annual, summer, and winter balance for each glacier and each year in our dataset.

Figure 3 .
Figure 3. Mass balance sensitivities (m w.e.year −1 ) to (a, b) temperature and (c, d) precipitation change for the 45 glaciers.In (b) and (d), the red line shows the mean sensitivity across the 45 glaciers, and the blue shading shows the standard deviation.

Figure 3 .
Figure 3. Mass balance sensitivities (m w.e.year −1 ) to (a,b) temperature and (c,d) precipitation change for the 45 glaciers.In (b) and (d), the red line shows the mean sensitivity across the 45 glaciers, and the blue shading shows the standard deviation.

Figure 4 .
Figure 4. Spatially heterogeneity in mass balance sensitivity of the 45 glaciers to (a,b) temperature change and (c,d) precipitation change.

Figure 4 .
Figure 4. Spatially heterogeneity in mass balance sensitivity of the 45 glaciers to (a,b) temperature change and (c,d) precipitation change.

Water 2019 , 21 Figure 5 .
Figure 5. Sensitivity of mass balance (m w.e.year −1 ) to (a) temperature and (b) precipitation change versus elevation for the selected subset of nine glaciers.

Figure 6
Figure 6 shows the ELA sensitivity for those nine selected glaciers.The continental glaciers (Muztag Ata and Siachen) have a relatively low ELA, relatively large accumulation zone, and relatively insensitive ELA to temperature changes.On the other hand, ELA of maritime glaciers are most sensitive to temperature change: a temperature change of +1 K can lift the ELA above the maximum glacier elevation, as shown for AX010, Parlung No. 390, and Baishui Glacier No. 1.A study based on ASTER data has also found that maritime glaciers in the Nyainqentanglha and Spiti Lahaul regions have lost a large portion of their accumulation areas during 2000-2016 [6].

Figure 5 .
Figure 5. Sensitivity of mass balance (m w.e.year −1 ) to (a) temperature and (b) precipitation change versus elevation for the selected subset of nine glaciers.

Figure 6
Figure 6 shows the ELA sensitivity for those nine selected glaciers.The continental glaciers (Muztag Ata and Siachen) have a relatively low ELA, relatively large accumulation zone, and relatively insensitive ELA to temperature changes.On the other hand, ELA of maritime glaciers are most sensitive to temperature change: a temperature change of +1 K can lift the ELA above the maximum glacier elevation, as shown for AX010, Parlung No. 390, and Baishui Glacier No. 1.A study based on ASTER data has also found that maritime glaciers in the Nyainqentanglha and Spiti Lahaul regions have lost a large portion of their accumulation areas during 2000-2016 [6].

Figure 6 .
Figure 6.Sensitivity of ELA to temperature change of +1 K.

Figure 6 .
Figure 6.Sensitivity of ELA to temperature change of +1 K.

Figure 7 .
Figure 7.The influence of parameter perturbations on the mass balance sensitivity to (a) +1 K temperature change and (b) +10% precipitation change.

Figure 7 .
Figure 7.The influence of parameter perturbations on the mass balance sensitivity to (a) +1 K temperature change and (b) +10% precipitation change.

2018, 10 , 21 Figure 8 .
Figure 8. Squared correlation coefficient (R 2 ) between the mass balance sensitivity and (a,b) annual precipitation, (c,d) degree day factor, and (e,f) CI, all calculated across the 45 glaciers and shown for a range of temperature change, ΔT, scenarios.

Figure 8 .
Figure 8. Squared correlation coefficient (R 2 ) between the mass balance sensitivity and (a,b) annual precipitation, (c,d) degree day factor, and (e,f) CI, all calculated across the 45 glaciers and shown for a range of temperature change, scenarios.

Water 2018 , 21 Figure 9 .
Figure 9. Spatial distribution of factors influencing the mass balance sensitivities: (a) degree day factors, (b) summer precipitation proportion (%) to annual precipitation (pink line the area of <60%), (c) change in the refreezing proportion (green circles indicate a refreezing decrease, while red circles indicate refreezing increase) under +1 K scenario, and (d) contribution of snow (accumulation), refreezing, and melt (ablation) to the mass balance sensitivity under +1 K temperature change.

Figure 9 .
Figure 9. Spatial distribution of factors influencing the mass balance sensitivities: (a) degree day factors, (b) summer precipitation proportion (%) to annual precipitation (pink line encloses the area of <60%), (c) change in the refreezing proportion (green circles indicate a refreezing decrease, while red circles indicate refreezing increase) under +1 K scenario, and (d) contribution of snow (accumulation), refreezing, and melt (ablation) to the mass balance sensitivity under +1 K temperature change.

Table 1 .
Initial range of parameter values used in the model optimization and statistics for each optimized parameter: min, max, mean and standard deviation across the 45 glaciers.
Note that σ is the standard deviation.

Table 2 .
Model sensitivity to parameter perturbations calculated as RMSE, R 2 and bias between the observed and modeled (parameter-perturbed) mass balance time series averaged across the 45 glaciers.The standard run is the simulation using all the initial (optimized) parameter values.