Spatial pattern of glacier mass balance sensitivity to atmospheric forcing in High Mountain Asia

Abstract

The complex topography and size of High Mountain Asia (HMA) result in large differences in glacier mass-balance variability and climate sensitivity. Current understanding of these sensitivities is limited by simplifications in past studies’ model structure. This study overcomes this limitation by using a mass-balance model to investigate the climatic mass-balance variability and climate sensitivity of 16 glaciers covering major mountain ranges in HMA. Generally, glaciers in the southeast have higher mass turnover while glaciers at the margins of HMA show higher interannual mass-balance variability. All glaciers are most sensitive to temperature perturbations in summer. The climatic mass balance of 15 glaciers is most sensitive to precipitation perturbations in summer or spring and summer, even if the seasonal accumulation peak is not in summer. Only one glacier's mass balance (Chhota Shigri Glacier) is most sensitive to precipitation perturbations in winter. Glaciers with high mass turnover and high summer-precipitation ratio are more sensitive to temperature perturbations. Sensitivity experiments reveal that besides the non-linearity of mass-balance temperature sensitivity, mass-balance precipitation sensitivity is non-linear as well. Furthermore, resolving the diurnal cycle of albedo, (re)freezing and the differentiation between liquid and solid precipitation are important to assess climate sensitivity of glaciers in HMA.

1. Introduction

The glaciers of High Mountain Asia (HMA) are important long-term hydrological reservoirs for municipal supply, irrigation, and hydropower (Pritchard, 2019; Immerzeel and others, 2020; Li and others, 2022). Furthermore, they pose a significant threat through glacial lake outburst floods (GLOFs) and glacier involved rock and ice avalanches (Carrivick and Tweed, 2016; Furian and others, 2021; Kropáček and others, 2021; Zheng and others, 2021).

Overall, HMA glaciers have been losing mass in recent decades (e.g. Neckel and others, 2014; Kang and others, 2015; Ke and others, 2017; Shean and others, 2020). Shean and others (2020) estimate a HMA-wide annual negative geodetic mass balance of (−0.19 ± 0.03) m w.e. a⁻¹ between 2000 and 2018. Rapid warming is enhancing glacier retreat (Yao and others, 2019). The warming rate of 0.34 K (decade)⁻¹ between 1961 and 2010 (Wei and Fang, 2013) is approximately twice as high as the global mean. For the period 2015 to 2100, further shrinkage between (29 ± 12) % (Representative Concentration Pathway 2.6) and (67 ± 10) % (Representative Concentration Pathway 8.5) is projected due to temperature and precipitation changes (Rounce and others, 2020). An ice volume loss of (−21 ± 1) % by 2100 is estimated even without further warming during the 21st century due to the imbalance of glacier mass balance in HMA (Miles and others, 2021).

Nevertheless, significant differences between the mass balances in different regions of HMA exist (e.g. Yao and others, 2012; Li and others, 2018; Bonekamp and others, 2019; Bhattacharya and others, 2021). For example, since 2000, glaciers in some regions had a nearly balanced or even positive mass budget. This phenomenon is called the Karakoram anomaly (Hewitt, 2005; Farinotti and others, 2020), which is also visible in parts of the Pamir and West Kunlun Shan (Gardelle and others, 2013; Brun and others, 2017; de Kok and others, 2020; Zhu and others, 2022). Reasons for the high variability of the average mass balance between regions are the influence of different large-scale circulation systems (e.g. Yao and others, 2012; Mölg and others, 2014; Fugger and others, 2022), local circulation systems (Xu and others, 2014; Curio and others, 2015; Ma and others, 2018), differences in snow to total precipitation ratios and precipitation seasonality (Bonekamp and others, 2019), and different glacier morphology like slope, aspect, topography and size (Yu and others, 2013; Brun and others, 2019).

Large-scale circulation systems which impact the glaciers in HMA are mainly the Indian Summer Monsoon (ISM), the westerlies and the East Asian monsoon system (Yao and others, 2012). Generally speaking, the ISM influences the spring- and summer-accumulation type glaciers (Benn and Owen, 1998; Fujita and Ageta, 2000; Yao and others, 2012). This influence decreases in the direction of the interior of HMA resulting in lower precipitation amounts and higher equilibrium line altitudes (ELAs). In the northwest (e.g. Pamir) precipitation rates are higher again (ELAs decrease) due to the influence of the westerlies (Yao and others, 2012).

Glaciers in wetter climates are more sensitive to temperature change (Oerlemans, 1997) and less sensitive to changes in precipitation (Yang and others, 2013). Ablation and temperature are non-linearly connected (Bolibar and others, 2022). With higher temperatures, melt rates increase, the melt season is prolonged and less precipitation falls as snow. In addition to the reduced mass input (snow), more shortwave radiation is absorbed due to lower albedo. Therefore, even more ablation is triggered (Oerlemans and Knap, 1998; Fujita, 2008b). Glaciers that have substantial amounts of accumulated snowfall in summer are extremely sensitive to the effect of less precipitation falling as snow (Qu and Hall, 2007; Fujita, 2008a; Ghatak and others, 2014; Johnson and Rupper, 2020; Arndt and others, 2021). The precipitation sensitivity in summer exists both for glaciers with high and low precipitation amounts. In the case of high precipitation amounts, the sensitivity to precipitation changes is controlled by the effect on accumulation. In the case of low precipitation amounts, it is controlled by the albedo influence on ablation (Fujita, 2008b). Glaciers in arid, cold climates are less sensitive to temperature changes (Ohmura and others, 1992; Braithwaite, 2008). Glaciers with large winter precipitation are less sensitive to changes in the liquid/solid fraction of precipitation (Kapnick and others, 2014). Such glaciers, as for example in the Pamirs, are more controlled by precipitation than by air temperature changes (Zhu and others, 2020).

Recent glacier mass balance sensitivity studies in HMA, confirm the higher sensitivity of summer accumulating glaciers (Sakai and others, 2015), a nonlinear temperature but linear precipitation sensitivity of glaciers (Wang and others, 2019), and that spatial differences in glacier response have a higher explanatory power than regional differences in climate change (Sakai and Fujita, 2017). The factors including summer temperature, temperature range, and summer to annual precipitation ratio are major indicators of spatial variability in mass-balance sensitivity (Sakai and Fujita, 2017).

These past studies, while informative, have been based on degree-day models (Wang and others, 2019), often at daily or coarser temporal resolutions, or focussing on the temperature sensitivity only (Sakai and Fujita, 2017). Instead, we apply the physically-based COupled Snowpack and Ice surface energy and mass balance model in PYthon (COSIPY, Sauter and others, 2020) to 16 different glaciers representing different climate regions and glacier types in HMA (see Fig. 1) with a consistent approach. Hourly meteorological variables of the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis ECMWF Re-Analysis fifth generation-Land (ERA5-L, Muñoz-Sabater and others, 2021) are statistically downscaled to the glaciers and used as forcing from January 2000 until September 2018. Total precipitation of ERA5-L is scaled in such a way that the simulated annual glacier-wide climatic mass balance fits the geodetic estimate of Shean and others (2020) for the individual glaciers. With this, we apply a physically uniform modelling chain, including the same climate forcing, mass-balance model and calibration dataset. We analyse the temporal and spatial variability and the sensitivity against seasonal and overall perturbations of temperature and precipitation.

Figure 1. Study region and location of the glaciers in High Mountain Asia with elevation represented in colours (Topographic-WMS, terrestris GmbH & Co, 2021) in the main map. Small inset maps with Randolph Glacier Inventory 6.0 (RGI Consortium, 2017) outlines of Abramov glacier (ABR), Batysh Sook glacier (SOO), Keli Yanghe source glacier (KYS), Urumqi Glacier No. 1 (UG1), Bayi Ice Cap (BIC), Muztagh Ata glaciers (MZA), Purogangri Ice Cap (PIC), Siachen glacier (SIA), Parlung No. 94 glacier (PL94), Chhota Shigri Glacier (CSG), Parlung No. 4 glacier (PL04), Guliya Ice Cap (GIC), Naimona'nyi glacier (NAI), Halji glacier (HAL), Yala Glacier (YAL) and Zhadang glacier (ZHA). RGI6 polygons for UG1, BIC, MZA and PIC have been merged. The RGI6 polygon of CSG is larger than the outline of CSG in other studies. However, we had to use the full RGI6 polygon because of the scaling to Shean and others (2020). In the case of UG1, the adjacent icefield is excluded because of an ice divide. Backdrops within the small inset maps are from Bing Maps (Microsoft, 2022).

The study region and the glaciers of interest are presented in the following Sec. 2. Section 3 describes COSIPY, the preprocessing of the climate and static data, and the executed simulations including the sensitivity studies. The results are presented in Sec. 4 and discussed in Sec. 5.

2. Glaciers of interest

This study includes glaciers in the Himalayas (southern and northern flank), the Nyainqentanglha, the Qilian Shan, the West Kunlun Shan, the Tian Shan, the Karakoram and the Pamir, representing all types of accumulation regimes proposed by Maussion and others (2014). We selected glaciers of different sizes because they are well-studied benchmark glaciers for specific regions. Shi and Liu (2000) classified the glaciers in HMA into three types: maritime (temperate), subcontinental (subpolar) and extreme continental (polar) glaciers mainly based on annual precipitation and annual and summer mean air temperatures at the ELA. All types are represented in this study. Furthermore, we wanted to include a glacier with a statistically-significant positive mass balance (see Fig. 4a in Shean and others, 2020).

Figure 1 displays the glaciers of interest including their outlines according to the Randolph Glacier Inventory 6.0 (RGI6, RGI Consortium, 2017). RGI6 glacier IDs, glacier area, maximum, minimum and mean elevation, geographic location and geodetic mass-balance by Shean and others (2020) for the period 2000 to 2018 are provided in Table 1. The bias-corrected (see Sec. 3.2) ERA5-L mean climate data of all glaciers are displayed in Table 2. Their long-term (2001–2018) mean annual cycles (except surface pressure) are displayed in Fig. S1, Fig. S2 and Fig. S3.

Table 1. Area-sorted properties of studied glaciers (RGI Consortium, 2017) and geodetic glacier mass balance according to Shean and others (2020) for the period 2000 to 2018.

Full glacier names are provided in the caption of Fig. 1.

Table 2. ERA5-L mean surface pressure p_sfc, air temperature at 2 m T ₂, relative humidity at 2 m RH ₂, incoming shortwave radiation Q_SWin, incoming longwave radiation Q_LWin, wind speed at 2 m U ₂, annual total precipitation TP and scaled annual TP.

Scaled TP is the precipitation required for the climatic mass balance to correspond to the geodetic mass balance of Shean and others (2020) for the individual glaciers. The values are the bias-corrected variables at the mean elevation of each glacier. The glaciers are sorted according to their area. Full glacier names are provided in the caption of Fig. 1. Uncertainties of the scaled TP amounts are based on the glacier mass balance uncertainties (see Table 1) reported by Shean and others (2020).

The glaciers are Batysh Sook (SOO, Kenzhebaev and others, 2017), Zhadang (ZHA, Huintjes and others, 2015), Urumqi Glacier No. 1 (UG1, Li and others, 2021), Yala (YAL, Stigter and others, 2018), Halji (HAL, Arndt and others, 2021), Bayi Ice Cap (BIC, Qing and others, 2018), Parlung No. 94 (PL94, Shaw and others, 2021), Keli Yanghe source (KYS), Naimona'nyi (NAI, Zhu and others, 2021), Parlung No. 4 (PL04, Shaw and others, 2021), Chhota Shigri (CSG, Azam and others, 2016), Abramov (ABR, Denzinger and others, 2021), Guliya Ice Cap (GIC, Thompson and others, 2018), Muztagh Ata (MZA, Holzer and others, 2015), Purogangri Ice Cap (PIC, Liu and others, 2019) and Siachen (SIA, Berthier and Brun, 2019). The literature cited to each glacier is a selection of recent publications only. To our knowledge, there has not been a specific mention or a study at the KYS. We selected this glacier because of its statistically-significant positive mass balance in West Kunlun Shan. According to Bing Maps (Microsoft, 2022) and Open Street Map (OpenStreetMap contributors, 2017), KYS is one of the source glaciers of Keli Yanghe, which in turn is a tributary of the Karakash river which drains into the Tarim Basin. There might be other local names for KYS we are not aware of.

Further, UG1 is explicitly named as an example for a subcontinental type glacier and GIC as an extreme continental type by Shi and Liu (2000). Eight glaciers have an area smaller than 4 km², 12 glaciers smaller than 22 km² and four glaciers are larger than 100 km² (RGI6 RGI Consortium, 2017). According to the Shuttle Radar Topography Mission (SRTM, Farr and others, 2007) digital elevation model (DEM), the glaciers are situated between 3900 m a.s.l. and 7500 m a.s.l. The hypsometry of the 16 glaciers in 50 m bands can be found in Fig. S4. Between 2000 and 2018, ZHA had the most negative mass balance with −1.06 m w.e. a⁻¹ and KYS the most positive with +0.43 m w.e. a⁻¹ (Shean and others, 2020). Three glaciers (KYS, MZA, SOO) had a positive mass budget, and one (SIA) was nearly balanced between 2000 and 2018.

3. Methods and data

In this chapter, COSIPY is introduced, the preprocessing of the climate forcing data is described and the setup for the performed simulations and analyses is explained.

3.1 COSIPY

The COupled Snowpack and Ice surface energy and mass balance model in PYthon (COSIPY) is a physically-based medium-complexity energy and mass-balance model (Sauter and others, 2020). The surface energy balance (SEB) is solved with an optimisation algorithm and coupled through the surface temperature T_s, surface melt M_sfc, and rain with an adaptive (non-equidistant) subsurface scheme. The modular design is intended to ensure maximum traceability, and easy optimisation and implementation of new parameterisations. The spatial and temporal resolution is freely scalable and dependent on the static and dynamic input data. The open-source model is coded in Python 3 and available on GitHub (https://github.com/cryotools/cosipy, last access: 28 November 2021) and in this study we apply COSIPY v1.3 (https://doi.org/10.5281/zenodo.3902191, last access: 28 November 2021). The SEB can be written as follows:

(1)$$Q_M = Q_{SWin} ( 1-\alpha ) + Q_{LWin} + Q_{LWout} + Q_H + Q_E + Q_G + Q_R$$

where Q_M is the available melt energy, Q_SWin is the incoming shortwave radiation, α is the albedo, Q_LWin is the incoming longwave radiation, Q_LWout is the outgoing longwave radiation, Q_H is the sensible heat flux, Q_E is the latent heat flux, Q_G is the glacier heat flux and Q_R is the sensible heat flux of rain. All terms on the right-hand side of Eq. (1) are input (Q_SWin and Q_LWin) or intermediate (α) variables, or parameterised (Q_LWout, Q_H, Q_E, Q_G and Q_R) using the optimisation algorithm. The residual of the optimisation algorithm results in Q_M only if T_s is at the melting point temperature.

COSIPY is a point model, i.e. does not consider ice-dynamics or any lateral exchange of energy and mass. Basal processes are not resolved. Therefore, the calculated mass-balance is defined as the climatic mass balance in accordance with Cogley and others (2011):

(2)$$b_{clim} = c_{sfc} + a_{sfc} + c_i + a_i = b_{sfc} + b_i$$

where c_sfc is surface accumulation, a_sfc is surface ablation, c_i is refreezing, a_i is subsurface melt, b_sfc is the surface mass balance and b_i is the internal mass balance. In COSIPY, c_sfc results from accumulated snowfall SF_c and in the case the surface temperature is <0 °C from deposition of water vapour from the atmosphere. Available melt energy from Eq. (1) and sublimation result in a_sfc. Refreezing of percolating water results in c_i, and subsurface melt M_sub resulting from penetrating radiation equals a_i.

All types of frozen precipitation equal SF_c because COSIPY does not consider direct sublimation during snowfall and processes associated with snowdrift. When total precipitation TP is used as input for COSIPY, instead of using snowfall directly, TP is separated into snow and rain as a function of air temperature at 2 m T ₂ using a logistic transfer function (Hantel and others, 2000). At 1 °C, the logistic transfer function results in 50 % snowfall and 50 % rain. The temperature range in which both phases occur is about 5 K (3.5 °C: 1 % snowfall, −1.5 °C: 99 % snowfall).

Energy fluxes are of positive (negative) algebraic sign towards (directed away from) the surface and presented in W m⁻². Mass-balance components are positive (negative) if they result in mass gain (loss) and presented in m w.e. Within COSIPY, the glacier-wide climatic mass balance B_clim is the sum of all point b_clim calculated by COSIPY at each glacier grid point (GGP). The number of GGPs depends on the resolution of the selected DEM. An advantage of the independence of the GGPs from each other is the easy implementation of the model on High-Performance Computing Clusters (HPCCs). The glacier area is thereby resolved by the GGPs derived from the DEM and the glacier outline. In addition, the annual glacier-wide climatic mass balance B_clim,a and the glacier-wide cumulative climatic mass balance B_clim,cum are used for the presentation of results.

In the manuscript the term refreezing refers to the refreezing of percolating meltwater and the freezing of percolating rain, subsurface melt and to a very small share the freezing of percolating condensation. The mass-balance year (MB-year) corresponds to the hydrological year (Cogley and others, 2011). It starts on 1 October and ends on 30 September. The year containing January in the mass-balance year (MB-year) determines its name. The mass turnover as used in this study is defined as the average of the absolute values of mean annual accumulation and mean annual ablation over the entire study period. In COSIPY, the albedo is parametrised according to Oerlemans and Knap (1998). Values for fresh snow, firn and ice albedo correspond to the value in Sauter and others (2020). We apply the parameters from Mölg and others (2012) for the parameters of time and height. For a comprehensive description of the model, please refer to Sauter and others (2020).

3.2 Atmospheric input data

ERA5-L is a subversion of ERA5 at higher resolution and used as model forcing (Muñoz-Sabater and others, 2021). From the forcing variables required for COSIPY (Table 3), only T ₂ and dewpoint temperature at 2 m T _d,2 are output variables of ERA5-L. The other forcing variables to COSIPY are also forcing variables for ERA5-L which implies that no substantial differences exist between ERA5-L and ERA5 data in these cases.

Table 3. COSIPY forcing variables, applied downscaling approaches to ERA5-L data and approaches to create the distributed fields (interpolation) on the glaciers.

A dash stands for no downscaling. The applied lapse rates are calculated for each glacier separately. The table is an adapted version of Table 2 of Arndt and others (2021).

First, the variables from the ERA5-L grid cell in which the centre point of the glacier is located are statistically downscaled to the elevation of the centre point. The height difference results from the modelled height of the ERA5-L grid cell and the mean elevation of the glacier. Afterwards, the variables are interpolated to the GGPs with the elevation of the applied DEM (see Sec. 3.3). The barometric formula is used for downscaling of surface pressure p_sfc to the mean elevation of the glacier as well as the interpolation to the GGPs. T ₂ and T _d,2 (intermediate variable to calculate relative humidity at 2 m RH ₂) lapse rates used for the downscaling are calculated from their all-time mean of the 121 ERA5-L grid cells around the glacier's centre point. The same lapse rate is used for the interpolation of T ₂ to all GGPs, while for RH ₂ no interpolation approach is applied. No downscaling is applied to incoming shortwave radiation Q_SWin and incoming longwave radiation Q_LWin. Also no interpolation is applied to Q_LWin. The radiation model after Wohlfahrt and others (2016) is used. It corrects Q_SWin considering slope and aspect for the interpolation of Q_SWin to the GGPs. The scaling of TP is used as calibration of the B_clim,a to the remote sensing derived geodetic glacier mass balance of Shean and others (2020). The derived wind speed at 2 m U ₂ of ERA5 and ERA5-L underestimates surface winds, especially in regions with high orography. One reason is the smoothing effect of valleys and mountains when using a 31 km grid (personal communication with ECMWF-Support, on 20 November 2018). Further, the surface wind from ERA5-L represents the grid cell's average wind (personal communication with ECMWF-Support, on 20 November 2018). Based on the experience from previous studies (Thiel and others, 2020; Arndt and others, 2021) we apply a constant scaling factor of two for all glaciers and no interpolation to the GGPs. Furthermore, no lapse rate is applied to TP resulting in the same amount of TP for each GGP. Nevertheless, SF_c differs at each GGP, since TP is partitioned into wet and solid precipitation depending on T ₂, to which, in turn, a lapse rate is applied.

3.3 Simulations and procedures

First, static files for all glaciers using a DEM and glacier outline are generated with the Geospatial Data Abstraction Library (Farr and others, 2007) using the COSIPY preprocessing utilities. The static files contain the binary glacier mask, the geographic location, the elevation, the aspect and the slope for each of the resulting GGPs. Outlines are derived from the RGI6 (RGI Consortium, 2017). The SRTM 1 arcsecond global product (Farr and others, 2007) is used as the DEM. There are ‘’no data” values for five glaciers (PL94, KYS, PL04, MZA, SIA). Hence, for these glaciers the post-processed ‘’Hole-Filled SRTM for the Globe Version 4” (Jarvis and others, 2008) is used. Next, ERA5-L hourly data (henceforth called forcing data) are downscaled (see Sec. 3.2) to the mean elevation of the glaciers. Simulations with coarser spatial resolution (fewer GGPs) revealed that a minimum of 150 GGPs leads to deviations < 0.1 m w.e. a⁻¹ compared to the highest possible spatial resolution. Only for the three largest glaciers (SIA, MZA, PIC) at least 300 GGPs were used. With fewer GGPs, the deviations increase and scatter (see Fig. 5, Arndt and others, 2021). With this in mind, we aggregated the static files for the best tradeoff between deviation and computational cost to resolutions between 60 m (SOO) and 1500 m (SIA, see Table S1).

For all following core COSIPY simulations (more than 1600), the HPCC of the Climate Geography lab, Humboldt-Universität zu Berlin, Germany was used. All simulations were executed using an hourly resolution for the period 1 January 2000 until 30 September 2018, while only the MB-years 2001 until 2018 are considered for all analyses and plots. The simulations were executed with a glacier-wide initial snow height of 0.2 m, a basal temperature of −1 °C, an initial top snow density of 300 kg m⁻³ and a bottom snow density of 500 kg m⁻³. The constants of the albedo and roughness parameterisation, and the snow transfer function are displayed in Table 4. The density of fresh snow is calculated depending on T ₂ and U ₂ in accordance with Vionnet and others (2012). For the initialisation of the density and temperature profiles, please refer to the values given by Sauter and others (2020). The resulting B_clim,a are compared with the geodetic mass balances of Shean and others (2020). If a glacier consists of multiple RGI6 outlines, results of Shean and others (2020) for the individual polygons have been weighted in relation to the total area of the glacier. Simulations are re-executed with scaling TP until the difference between both mass budgets is ≤0.03 m w.e. a⁻¹. Between 4 and 10 simulations per glacier had to be performed to find the single scaling factor for each glacier. Resulting scaling factors vary between 0.48 and 1.59 (see Table S1). The calibration procedure compensates for all uncertainties in climate forcing data, as well as model parameters and any processes not covered by the model (see Sec. 5.1.1). A prerequisite of this study is to use the same forcing data with the same approach for bias correction, the same model setup and the same calibration data for all glaciers throughout to exclude any differences that may arise from different data, approaches or models.

Table 4. Varied climate forcing variables and model parameters of the sensitivity experiments.

To quantify the climate sensitivity of the glaciers, we calculated seasonal sensitivity characteristics SC_sea following Oerlemans and Reichert (2000) and uniform sensitivity characteristic SC_uni similar to Schuler and others (2005). For the SC_sea, we first applied T ₂ perturbations for each glacier to obtain an equilibrium B_clim,a ≤ |0.04| m w.e. for the whole study period. The required perturbations are within the recommended range of ±2 K (Oerlemans and Reichert, 2000) with values between −1.32 K (PL04) and 1.07 K (MZA, see Table S1). Subsequently, 24 simulations with single monthly temperature perturbations (each calendar month 0.5 K and −0.5 K) and 24 simulations with single monthly precipitation perturbations (each calendar month 10 % and −10 %) are executed. The 12 temperature pairs and 12 precipitation pairs are the indexes which are plotted for each glacier and display the deviation of B_clim,a to the monthly perturbations.

For the SC_uni, we executed uniform temperature perturbations between −1.5 K and 2.5 K in 0.5 K steps and uniform precipitation perturbations between −50 % and 50 % in 10 % steps. Unlike the SC_sea, the reference runs for the SC_uni are the simulations calibrated to Shean and others (2020).

4. Results

4.1 Simulated energy and mass balance

The simulated annual glacier-wide climatic mass balance B_clim,a, glacier-wide cumulative climatic mass balance B_clim,cum, annual accumulated snowfall SF_c, annual rain and runoff of the 16 studied glaciers for the MB-years 2001–2018 are displayed in Fig. 2. The MZA and KYS have only positive B_clim,a, PIC, PL04, PL94, ZHA, GIC have only negative B_clim,a and the rest have positive and negative B_clim,a. The glaciers with the most negative B_clim,cum (except ZHA) are glaciers that have positive and negative B_clim,a, in contrast to GIC and PIC on the Tibetan Plateau which have only negative B_clim,a but have not the most negative B_clim,cum over the whole study period. B_clim,a of GIC, PIC, ZHA and KYS – all located centrally within HMA – have the smallest interannual variability according to the coefficient of variation c_v, with values < 0.15 (see Fig. 2). The detrended standard deviation of B_clim,a depends on the amount of B_clim,a. Therefore, c_v is defined as the ratio between the detrended standard deviation of B_clim,a and the mean of the annual values of the mass turnover similar to Maussion and others (2014). B_clim,a of glaciers in the Tian Shan occur with a high interannual variability with c_v > 0.40. B_clim,a of the glaciers in the Himalayas are also highly variable with a mean c_v > 0.35 including the highest interannual variability of NAI. B_clim,a of SIA also has a high interannual variability (c_v = 0.44), while the B_clim,a glaciers in the Pamir, Pamir-Alay and Qilian Shan have a medium interannual variability (c_v = 0.28,  0.29,  0.31). B_clim,a of Parlung glaciers have quite different interannual variability with c_v values of 0.25 and 0.39.

Figure 2. Annual glacier-wide climatic mass balance B_clim,a (black), accumulated snowfall (blue), rain (light blue) and runoff (red) of mass-balance years 2001–2018 (left y-axis); glacier-wide cumulative climatic mass balance (black line, right y-axis). Standard deviation σ of B_clim,a and coefficient of variation c_v (σ of B_clim,a divided by mass turnover). Full glacier names are provided in the caption of Fig. 1.

In general, the glaciers at the southern and eastern margin have the highest mass turnover and the most negative B_clim,cum, while the ZHA with the most negative B_clim,cum is not located at the margin. The ABR and UG1 contradict this pattern. Located at the northwestern and northern margin respectively, these glaciers have a high mass turnover as well.

The mean monthly (2001–2018) values of SF_c, rain and runoff are displayed in Fig. 3. Except for KYS, all glaciers have their ablation peak in August (12 glaciers) or July (3 glaciers). The ablation is very low at KYS with runoff amounts in the peak summer months of ~0.01 m w.e. month⁻¹ and a mean value (MB-year 2001–2018) of only 0.007 m w.e. month⁻¹ in July. Using SF_c as criterion for the accumulation-season classification, the glaciers PL04, PL94, SOO, UG1 would be spring-accumulation types. The glaciers ABR, SIA and CSG are hybrid winter/spring-accumulation type glaciers. The glaciers KYS, GIC, MZA, NAI, PIC, BIC, and YAL are summer-accumulation type glaciers. ZHA is a spring to summer-accumulation type glacier. HAL accumulates solid precipitation in winter, spring and summer and only in November and December the amount is close to zero.

Figure 3. Mean monthly (mass-balance years 2001-2018) accumulated snowfall (blue), rain (light blue) and runoff (red). The sum of percolating surface melt, subsurface melt and rain reaching the snow/ice interface is named ‘’runoff”. Full glacier names are provided in the caption of Fig. 1.

The glacier-wide overall mean energy fluxes are presented in Fig. 4. Net shortwave radiation Q_SWnet is the greatest energy source for all glaciers. GIC, and ZHA located at the interior of HMA receive the highest Q_SWnet compared to the other glaciers. The glaciers with the most positive B_clim,a KYS receive the lowest energy input from Q_SWnet because of high glacier-wide albedo. Furthermore, glaciers in the same mountain ranges receive energy from Q_SWnet in the same range. Glaciers in the Himalayas receive between 57 W m⁻² and 64 W m⁻², the Parlung glaciers between 47 W m⁻² and 54 W m⁻², and glaciers in the Tian Shan, Pamir and Pamir-Alay between 44 W m⁻² and 51 W m⁻². Net longwave radiation Q_LWnet is the most important energy sink for all glaciers. The glaciers with the highest energy loss compared to all other glaciers are GIC, PIC and ZHA located at the interior of HMA. The glacier-wide overall mean Q_E of each glacier is an energy sink, while Q_H is an energy source for some glaciers and an energy sink for the others. The sign of Q_E and Q_H is determined by the temperature gradient and the gradient of the mixing ratio of water vapour, respectively. Because of the transport of latent energy through surface melt and rain into the snow and resulting latent heat release through refreezing and the energy transfer via penetrating shortwave radiation in the uppermost layers, Q_G is always positive and directed to the surface in the overall mean. Q_R plays a minor role and is certainly always positive. Q_SWnet, Q_LWnet and the closing term of the energy balance Q_M are the dominant energy components of all glaciers.

Figure 4. Overall mean net shortwave radiation (yellow), net longwave radiation (blue), sensible heat flux (dark green), latent heat flux (light green), glacier heat flux (black), sensible heat flux of rain (orange) and available melt energy (red). Full glacier names are provided in the caption of Fig. 1.

Figure 5 presents the B_clim,a components. The most important mass flux is M_sfc for all glaciers with the exemption of KYS, which has the least negative annual rate with −0.44 m w.e. a⁻¹. The most negative mean annual value is YAL with −2.83 m w.e. a⁻¹. With the exemption of MZA, refreezing is higher than SF_c, while it has to be kept in mind that refreezing depends on the available percolation water through M_sfc, rain, subsurface melt and to a small share condensation of water vapour at the surface. Approximately 85 % of percolating water freezes in the case of the KYS within the snowpack (see Table S2). YAL has the lowest ratio with 29 %. Most glaciers (11 out of 16) show values between 29 % and 43 %. Besides the exceptionally high ratio of KYS, the four glaciers MZA (63 %), SIA (59 %), SOO (58 %) and NAI (51 %) have above-average ratios. With the exception of KYS and SOO these glaciers are the glaciers which extend over an elevation range > 1500 m. YAL – located on the southern flank of the Himalayas – has the highest amount of SF_c. M_sub and sublimation have a minor influence on the mass budget. Despite high Q_E of some glaciers, sublimation plays a minor role in the mass balance due to the high energy demand to sublimate a certain amount of snow or ice in comparison to the demand for melting the same amount of snow (latent heat of sublimation: 2.832 10⁶ J kg ⁻¹ K ⁻¹ $\gg$ latent heat of fusion: 3.337 10⁵ J kg ⁻¹ K ⁻¹).

Figure 5. Overall mean accumulated snowfall (blue), rain (light blue), surface melt (red), refreezing (cyan), subsurface melt (orange), sublimation (grey) and annual glacier-wide climatic mass balance (black). Full glacier names are provided in the caption of Fig. 1.

The b_clim,a profiles in 50 m bands for the MB-years 2001–2018 are displayed in Fig. S5. No 50 m band has a higher annual accumulation rate than 1.5 m w.e. a⁻¹, while the bands with the highest ablation rates reach 6 m w.e. a⁻¹. The SOO has the lowest ELA with 4106 m a.s.l. and NAI the highest with 6002 m a.s.l. A clear north/south ELA pattern is visible. All three glaciers north of 39.5 $^\circ$N have ELAs below 4300 m and all five glaciers south of 30.5 $^\circ$N ELAs above 5200 m.

To identify geographical patterns, we tested Spearman's rank correlation of B_clim,a, SF_c, runoff and mass turnover with latitude, longitude and the composite of latitude and longitude (longitude + [90-latitude]) which describes the northwest/southeast component of HMA. Figure S6 supports the pattern with a more negative B_clim,a and higher mass turnover on the diagonal from northwest to southeast at the significance level 0.01.

4.2 Mass balance sensitivity

The calculated 24 indexes of the seasonal sensitivity characteristics SC_sea are displayed in Fig. 6. All glaciers have their seasonal peak sensitivity to monthly temperature perturbations in summer, most of them in July (BIC, PL94 in August and GIC in June). Most of the glaciers are also most sensitive to changes in precipitation amount in summer, but the precipitation sensitivity is less concentrated to summer months. They also react to precipitation perturbations in late winter and spring. ZHA with the most negative B_clim,a is also the most sensitive glacier to both temperature and precipitation changes in a single month.

Figure 6. Seasonal sensitivity characteristics (SSC) after Oerlemans and Reichert (2000). Red bars are the dependence of the annual glacier-wide climatic mass balance B_clim,a on monthly temperature perturbations of 1 K and blue the dependence of B_clim,a on monthly total precipitation perturbation of 10 %. Full glacier names are provided in the caption of Fig. 1.

The specific pattern of the glaciers, which extend over a high elevation range and have a low-elevated glacier tongue related to their elevation range, is visible in the SC_sea as well. SIA, MZA, PL04, NAI and CSG are the least sensitive to monthly temperature or precipitation perturbations. In this context, it is worth comparing the two pairs of nearby glaciers. HAL (elevation range: 400 m) and PL94 (elevation range: 577 m) are much more sensitive to temperature changes than NAI (elevation range: 1686 m) and PL04 (elevation range: 1293 m). BIC and PIC react similarly, keeping in mind that they are very different in size. KYS is most sensitive to changes in temperature and precipitation perturbations in July. For no other glacier is the difference in sensitivity between the most sensitive month and the second most sensitive month so pronounced.

The annual cycle of the seasonal sensitivity to precipitation perturbations is similar to the long-term (2001–2018) annual cycle of SF_c (see Fig. 3). Some differences are evident, though, such as SF_c in July and August for SIA is lower than in January to March, while the sensitivity to precipitation perturbations is highest in summer. Moreover, the precipitation sensitivity of the spring-accumulation type glaciers (PL94, PL04, and hybrid spring/summer accumulation type ZHA) is shifted to summer. Noteworthy is the sensitivity to precipitation perturbations in February and March of GIC, while the long-term mean of SF_c is very low in these months. PIC, BIC and UG1 show a peak of SF_c and of precipitation SC_sea in June. CSG is the only glacier that reacts most sensitively to precipitation perturbations in winter. The peak of the precipitation sensitivity and also accumulation by SF_c is in February.

Similar to Schuler and others (2005) the SC_uni are shown in Fig. 7. Noteworthy is the shape of the SC_uni temperature curves of YAL, ZHA and BIC. It seems that they are currently in a temperature regime, which is the most sensitive for them, i.e. the slope of the curve is the highest at 0 K perturbation. Moreover, the three glaciers are the most sensitive glaciers to a uniform perturbation of 2.5 K and among the five most sensitive glaciers to a uniform TP reduction of 50 %. According to the slope of the lines, the glaciers SIA and MZA, which extend over the highest elevation range, are the least sensitive to uniform temperature perturbations, while NAI which also extends over a large elevation range is one of the most sensitive glaciers. The reasons are possibly the different climate background (higher summer T ₂, higher summer precipitation ratio) of NAI in comparison to SIA and MZA (see Sec. 5.3.2). SIA and MZA react more similarly to the other glaciers in terms of uniform TP perturbation than T ₂ perturbations. Noteworthy is the B_clim,a response of nearly −8 m w.e. a⁻¹ of YAL to a TP perturbation of −50 %. The patterns are similar when the starting point for the perturbations are the temperature adjusted simulations of the SC_sea to reach a zero B_clim,cum (see Fig. S7) and not the ones scaled to Shean and others (2020).

Figure 7. Uniform sensitivity characteristics SC_uni. Displayed are the annual glacier-wide climatic mass balance B_clim,a changes to the simulations without perturbations. Full glacier names are provided in the caption of Fig. 1.

5. Discussion

5.1 Uncertainties and validation

5.1.1 Forcing data, model uncertainties and simplifications

Considerable uncertainty is related to the reanalysis forcing. We used linear lapse rates for T ₂ and did not correct the near-surface temperature due to the influence of katabatic winds, as parameterisations for their representation in models are still under development (e.g. Shaw and others, 2021). The error in Q_SWin due to the height difference between the ERA5-L grid cell and the glacier centre due to atmospheric attenuation is estimated to be small. In this study the ERA5-L Q_LWin output is used directly, but better modelling of cloud cover and Q_LWin is an important objective (Kok and others, 2019). Scaling of U ₂ was based on prior studies, but introduces additional uncertainty into the computations for newly analysed glaciers.

The forcing variable TP is used as a spatially-uniform calibration variable. TP is the most critical forcing variable. Without local measurements of TP it is not possible to use TP directly (e.g. Immerzeel and others, 2015; Gao and others, 2018; Hamm and others, 2020; Arndt and others, 2021). To account for orographic effects on precipitation, high resolution precipitation modelling (Hamm and others, 2020) or at least information about region-specific lapse rates (Machguth and others, 2009) would be needed. A larger amount of SF_c is simulated in the higher-elevated areas compared to the lower-elevated areas of the glaciers due to the lapse rate of T ₂ and the fact that T ₂ separates TP into either rain or SF_c.

Of the COSIPY subroutines, B_clim,a is most sensitive to the albedo parameterisation (e.g. Arndt and others, 2021), as is typical for such models (e.g. Zhu and others, 2018). Different climate regimes and SF_c and T ₂ patterns lead to different constants for the albedo parameterisation (Oerlemans and Knap, 1998) for individual glaciers. Nevertheless, due to limited data it is not possible to use glacier-specific parameters. B_clim,a is sensitive to the separation of TP into solid and liquid based on the logistic transfer function, which is highly dependent on T ₂ and the parameter ‘’centre snow transfer” (50 % snow and 50 % rain).

Blowing and drifting snow, debris cover and mass accumulation by avalanches from surrounding mountain walls are not considered but can be important processes for mass-balance modelling. Blowing and drifting snow lead to mass input or output and impact albedo. Further, it changes the water vapour gradient near the surface and therefore ignoring these processes may overestimate sublimation due to missing saturation from drifting and blowing snow (Bintanja and Reijmer, 2001; Zhu and others, 2018). Debris-cover enhances melting (lower albedo) until a certain thickness and reduces melting (isolating effect) above the threshold (Östrem, 1959). Further, it changes the energy balance because surface temperatures above the melting point temperature are possible (Fugger and others, 2022). SIA is covered between 3 % and 13 % with debris (Agarwal and others, 2017; Kumar and others, 2020). It is assumed to be the glacier with the highest relative debris cover in this study. The contribution of avalanches can be especially important for glaciers with accumulation areas surrounded by steep valley walls (e.g. Laha and others, 2017) and a winter accumulation peak. Nevertheless, due to limited data it is not possible to use glacier-specific parameters within the study.

Due to the lack of consideration of ice dynamics and inconsistent information on area changes, no changes of hypsometry and area could be considered, although calculated sensitivities change as glacier geometry changes. However, since the simulations do not cover multiple decades such impacts are limited.

5.1.2 Model uncertainty experiment

The package Statistical Parameter Optimization Tool for Python (SPOTPY, Houska and others, 2015) was used to quantify the global sensitivity of B_clim,a of the 16 glaciers based on the Fourier amplitude sensitivity test (FAST) of Saltelli and others (1999). We divided the experiment into the uncertainties of climatic forcing data and model parameters, using the ranges as shown in Table 4. The spatial resolution was decreased so that the glaciers were represented by only 32 (BIC) to 80 (SIA) GGPs because of computational cost. This resulted in a maximum B_clim,a deviation of −0.11 m w.e. a⁻¹ (PIC and GIC) in comparison to the reference run. Furthermore, only 3 years (2001–2003) were simulated for each glacier. Tests showed that the surface roughness parameters have a small influence. Therefore, they were excluded from further experiments. One thousand simulations per glacier were carried out with varying climate forcing data (7 variables) and 700 simulations with varying model parameters (6 parameters).

The frequency distributions of the forcing data and model parameter experiments are shown in Figs. 8 and 9. Standard deviations vary between 0.38 m w.e. a⁻¹ and 1.94 m w.e. a⁻¹ for climate forcing variables and between 0.46 m w.e. a⁻¹ and 2.29 m w.e. a⁻¹ for model parameters. Zhu and others (2018) calculated a change in mass balance of 0.92 m w.e. a⁻¹ for the increase of one single albedo parameter (+10 % albedo fresh snow) at ZHA. The results here further confirm the results of Zhu and others (2018) regarding the extremely high sensitivity of ZHA compared to PL94. Wang and others (2019) calculated a mass balance response (averaged over 45 glaciers) of −1.59 m w.e. a⁻¹ from varying temperature lapse rates by −20 % and the least response of −0.45 m w.e. a⁻¹ by varying (−20 %) the snow/rain threshold in their parameterisations. Within both studies, these were only changes due to variation of a single parameter; in the case of Wang and others (2019) with monthly resolution and in the case of Zhu and others (2018) with daily resolution as opposed to the hourly resolution applied here. The resolution of the diurnal cycle has an important influence on model results, especially affecting precipitation phase partitioning, refreezing and albedo parameterisation. Most of the frequency distributions of the sensitivity experiments presented here are not normally distributed. As soon as B_clim,a becomes positive, the sensitivity to most parameters decreases drastically.

Figure 8. Frequency distributions of the forcing data sensitivity experiment. The x-axis shows the annual glacier-wide climatic mass balance of the 1000 simulations per glacier in m w.e. a⁻¹. The red line displays the reference simulation and σ is the standard deviation of the simulations. The y-axis displays the frequency. Please note the different scaling of the y-axes. Bins are also automatically scaled per glacier and not uniform. Full glacier names are provided in the caption of Fig. 1.

Figure 9. Frequency distributions of the model parameter sensitivity experiment. The x-axis shows the annual glacier-wide climatic mass balance of the 700 simulation per glacier in m w.e. a⁻¹. The red line displays the reference simulation and σ is the standard deviation of the simulations. The y-axis displays the frequency. Please note the different scaling of the y-axes. Bins are also automatically scaled per glacier and not uniform. Full glacier names are provided in the caption of Fig. 1.

Despite the high uncertainty of the equifinality solution between forcing data and model parameters, we can analyse and rank geographic differences of energy and mass balance components and mass balance sensitivity of the presented glaciers in HMA. However, absolute numbers should be evaluated with caution. The results should be treated as an estimation (order of magnitude) for the energy and mass balance patterns. In future studies, using additional independent datasets for calibration of B_clim,a could reduce the resulting equifinality problem (Barandun and others, 2018).

5.1.3 Validation

To ensure reliability of the results, we aim to evaluate the downscaled forcing variables, energy fluxes and resulting glacier mass balance. However, the study suffers from inconsistency of in situ observations between sites. Consequently, we demonstrate that the forcing and results are reasonable according to independent datasets.

We applied a scaling factor of 2 to the ERA5-L U ₂ because U ₂ is underestimated in the dataset. A comparison between the mean, scaled U ₂ and the mean measured U ₂ at PL94, ZHA, CSG and YAL (Table S3) indicates that the values are in the right order of magnitude.

The glaciers with the lowest and highest scaled annual TP are GIC and YAL, respectively. For GIC, Muhammad and Tian (2020) measured a TP amount of 185 mm a⁻¹ for one year compared to the 2000–2018 mean of 246 mm a⁻¹ here. The scaled TP amount of YAL is 2079 mm a⁻¹, which is higher than reported by nearby in situ observations (cf. Stumm and others, 2021). However, taking into account the possible undercatch of snow (Rasmussen and others, 2012), the average value of 1867 mm a⁻¹ (Sharma and others, 2020, data: Department of Hydrology and Meteorology, Kathmandu) for the whole of Nepal, and the fact that YAL is located on the southern slope of the Himalayas, this does not seem to be unrealistic. The general geographical pattern of scaled TP of the glaciers in this study is in accordance with the spatial precipitation variability in Fig. 5 of Maussion and others (2014).

The mean energy fluxes Q_SWnet, Q_LWnet, Q_H and Q_E for ZHA are all between the maximum and minimum values given in four other studies summarised in Table 7 in Zhu and others (2018). Mölg and others (2012) and Huintjes (2014) include shortwave penetrating radiation in the calculation of Q_G, in contrast to this study. Adjusting Q_G to include shortwave penetrating radiation in this study results in values for Q_G being similar to both studies and with the same sign. The resulting Q_M value of 21.9 W m⁻² for ZHA is close to the mean value (19.4 W m⁻²) of the four other studies (maximum: 30.6 W m⁻², minimum: 13.7 W m⁻²).

The interannual variability of the mass balance series is validated in Fig. S8, Fig. S9 and Fig. S10 with data of the World Glacier Monitoring Service (WGMS, 2021). Mass balance measurements exist for eight of the 16 glaciers, but only from three glaciers does a time series of more than 13 years exist. Therefore, only 3 glaciers could be statistically compared. For two of the three glaciers the variance of the interannual variability can be statistically significantly (p-value < 0.01) well represented (r ² = 0.59 and r ² = 0.4) by the modelling chain. For PL94, the result is not significant with r ² = 0.2. These results are expected, since no local-specific bias corrections or forcing data were used. The test shows that the general annual pattern of the different climate regions can be reproduced. Simulated and measured mass balances of the five glaciers with short times series are compared as B_clim,a (Fig. S11) and averaged mass balance profiles (Fig. S12). While the former show large differences for some glaciers, the latter show that their altitudinal patterns are reasonably similar even though no sophisticated parameter and bias-correction inference strategy (e.g. Rounce and others, 2020) with multiple indicators (e.g. Barandun and others, 2018) was applied.

Figure S13a shows the comparison of the mean B_clim,a simulated over the entire study period with the estimated elevation changes of Hugonnet and others (2021). It shows a good (r ² = 0.68) agreement of the change in the different mountain ranges. However, the dataset of Hugonnet and others (2021) is not strictly independent, since both the estimate of Shean and others (2020) and Hugonnet and others (2021) are derived from ASTER DEMs. Furthermore, we compared simulated ELAs and ELAs derived by Miles and others (2021) of the studied glaciers in Fig S13b. Again, the significant correlation (r ² = 0.79) shows the robustness of the modelling approach.

5.1.4 Uncertainty related to refreezing

Overall, we simulate a high ratio of refreezing to percolating water (surface melt + rain + subsurface melt + condensation). Fujita and Ageta (2000) modelled a rate of 22 % at the Xiao Dongkemadi. Xiao Dongkemadi glacier is most similar to the PIC in terms of its characteristics. For PIC we obtain a value of 41 %. Li and others (2018) simulated a rate of 25 % for Qiangtang No.1 Glacier, which is also most similar to the PIC. Yang and others (2013) calculate a refreezing percentage (without rain) of 9 % for PL94 compared to 32 % in this study. All three studies use daily values. In addition, these studies are point-scale studies, which hamper the comparison with the glacier-wide values in this study. Veldhuijsen and others (2021) investigated the meltwater refreezing rate of snow in the Langtang catchment, Himalayas. The refreezing decreases by 84 % when moving from hourly to daily resolution according to their findings. If a reduction in refreezing by 84 % (ratio decreased from 0.21 to 0.05) between hourly and daily forcing data was hypothetically applied to this study, the refreezing rates of PIC and PL94 would be reduced to 6 % and 5 %, respectively, which would result in refreezing well below the values based on daily data provided by Fujita and Ageta (2000), Li and others (2018) and Yang and others (2013).

The diurnal melt cycle is important, especially when the surface temperature is at the melting point temperature only for a short period each day, while deeper snow-layers are below the melting point temperature throughout. In such cases, snow-layers provide cold content for refreezing, which cannot be resolved based on daily data because no melt would occur in such cases. Higher melt rates based on hourly forcing are partly compensated by higher refreezing rates. Therefore, the difference in B_clim,a using daily forcing would probably not result in similar amounts as the differences in refreezing between hourly and daily resolution.

To further investigate the refreezing component we determined the long-term mean monthly mass fluxes of three glaciers (Fig. S14) and their annual portion of refreezing (Fig. S15). Refreezing occurs only when surface melt, rain or both are available. Furthermore, the spatial plot of KYS reveals that there is a spatial gradient of refreezing in the accumulation area from lower values at the highest locations to higher values at lower elevation. In the highest parts the refreezing is low because the surface melt is low. Only at two GGPs is there a low refreezing amount. These grid points are the only ones that are snow-free for longer periods during summer. In the case of NAI there is refreezing from the start of the ablation season until its end. NAI has a large, perennially snow-covered accumulation area. In this high elevation area the ratio of surface melt that refreezes within the snow layers (also in layers from previous years) is high. If rain occurs, much of the percolation rainwater also freezes. In the case of HAL the seasonal refreezing peaks in June while the ablation peak is in August. Diurnal plots of important energy- and mass-balance fluxes regarding refreezing at KYS and NAI at specific GGPs are shown in Fig. S16 and Fig. S17 to illustrate that the model produces physically consistent output. In conclusion, the refreezing rates in this study are high, but they are within a realistic range and can be well justified.

KYS stands out because barely any runoff is simulated resulting in the very large refreezing to percolating water rate of 85 %. KYS was selected because it is an example glacier within the West Kunlun Shan with one of the the most positive mass balances according to Shean and others (2020, see Fig. 4a). Nevertheless, we cannot evaluate the results of a glacier with such a positive mass budget between 2001–2018 because of a lack of other studies on such positive mass balances in HMA. Therefore, the results of KYS should be taken as an estimate for a glacier with such a positive mass balance.

5.2 Variability of energy and mass balance

Net shortwave radiation is typically the main energy source for mountain glaciers located at low to medium latitude (e.g. Greuell and Smeets, 2001; Oerlemans and Klok, 2002). Further, the expectation of higher sublimation rates for more continental type glaciers (e.g. Zhu and others, 2015) is reflected in the results, while the rate decreases again towards the western margin of HMA. With the percolation of surface water through the snow layers, latent energy is transported and released by refreezing. In addition, penetrating shortwave radiation heats up the uppermost layers, and, together with the latent energy from refreezing, leads to the positive glacier-average Q_G. A northwest–southeast pattern is most evident in B_clim,a and total mass turnover (see Fig. S6).

The definition of the accumulation season of the glaciers depends on whether TP or SF_c is taken as a classifying variable. Some glaciers are classified differently here than in other studies since we refer to SF_c only. Considering TP, most of the glaciers (except SIA, CSG, PL04, PL94) are summer-accumulation type glaciers, but the pattern changes considering only SF_c. For example, HAL, ZHA and PIC are not clear summer-accumulation type glaciers using SF_c as a classifying variable. PL04 and PL94 would be more spring than spring/summer accumulation type glaciers, which also better fits the results of Yang and others (2013) and Jouberton and others (2022). The percentage of TP occurring as SF_c in this study fits the pattern of Maussion and others (2014), while the datasets are fully independent.

5.3 Sensitivity

5.3.1 Seasonal sensitivity characteristics

Within this study, ABR was the only glacier studied by Oerlemans and Reichert (2000). The results point towards higher overall T ₂ sensitivity of ABR compared to Oerlemans and Reichert (2000). Different applied models and climate forcing are reasons for the different ranges. In the following, the focus is on seasonal patterns and the differences between studied glaciers, while absolute values of sensitivity are not further interpreted.

In comparison with the glaciers studied by Oerlemans and Reichert (2000), T ₂ sensitivity is more restricted to summer months. Only CSG has a distinct winter maximum in the sensitivity to TP. Sensitivity is largely restricted to the summer season for both T ₂ and TP perturbations for the three most continental type glaciers PIC, KYS and GIC and in accordance with similar findings for continental-type glaciers by Oerlemans and Reichert (2000). Further, sensitivity increases when the climate is wetter (Oerlemans, 1997; Fujita, 2008a) as a general pattern. Moreover, SC_sea reveals that B_clim,a is more sensitive to TP perturbations in summer than in winter, which reflects the strong influence on ablation through albedo. MZA, HAL and PL94 have a peak in summer TP sensitivity, while no clear peak in the monthly amount of SF_c in summer is present. All glaciers with a high summer precipitation amount are very sensitive to summer changes in T ₂, which is in accordance with Naito (2011). Total mass turnover is important for the SC_sea as well, which is supported by the statistically significant correlation (r ² = 0.38, p-value = 0.01) between the mass turnover and overall SC_sea (Fig. S18). Temperature increase in July would have a great influence on KYS in which case the ablation would increase tremendously. In months other than July, a temperature perturbation of 0.5 K is not sufficient to lead to negative B_clim,cum at KYS. Only in August and to some extent June, T ₂ temperature perturbations have an influence on B_clim,a.

5.3.2 Uniform sensitivity characteristics

Figure 7 reveals that neither the sensitivity to T ₂ nor to TP is linear (cf. Oerlemans and Fortuin, 1992; Fujita, 2008a). This partially contradicts Wang and others (2019), who found a linear relationship between TP perturbations and mass balance response. The combination of hourly forcing and the resolved albedo variability are probable reasons for the non-linear TP response.

Topographic and climatic characteristics and different sensitivity indexes are displayed in Table 5. Sakai and Fujita (2017) identified summer temperature, temperature range and summer precipitation ratio as the variables with the most explanatory power for the B_clim,a sensitivity to T ₂. The B_clim,a response of glaciers with a summer TP ratio < 50 % (SOO, PL94, PL04, CSG, ABR, MZA, SIA) results in a low T ₂ response. ZHA, UG1, YAL and HAL have summer precipitation ratios > 50 %. Nevertheless, from the four glaciers UG1 and HAL have high ΔT ₂ ranges resulting in a low T ₂ response, while ZHA and YAL have low ΔT ₂ ranges resulting in a high T ₂ response < −2.95 m w.e. a⁻¹ K⁻¹.

Table 5. Topographic and climatic characteristics and sensitivity indices.

Elevation range (ER); summer temperature (T ₂ JJA) calculated from the long-term (MB-year 2001–2018) mean of June, July and August; temperature range (ΔT ₂) calculated from the difference of the long-term mean of the warmest and the coldest month; summer total precipitation ratio (TP ratio) calculated as the percentage of June, July, August precipitation to annual precipitation; annual glacier-wide climatic mass balance sensitivity to a +1 K T ₂ perturbation (ΔB _T2) derived from uniform sensitivity characteristic SC _uni; annual glacier-wide climatic mass balance sensitivity to a +20 % TP perturbation (ΔB _TP) derived from SC_uni. The glaciers are sorted according to their area. Full glacier names are provided in the caption of Fig. 1.

GIC and PIC have summer T ₂ of −5 °C and −0.9 °C, respectively, resulting in low and medium T ₂ sensitivity. BIC has the highest T ₂ response with a value of −3.4 m w.e. a⁻¹ K⁻¹. Only high summer T ₂ and summer precipitation ratio > 50 % in combination can explain this high sensitivity to some extent. KYS has a low summer T ₂ resulting in the very low sensitivity even though the summer TP ratio is 60 %. The influence of summer TP ratio on T ₂ sensitivity is supported by an r ² of 0.29 (p-value = 0.03). A multiple regression model between the three explanatory variables summer T ₂, ΔT ₂ range and summer TP ratio and the T ₂ sensitivity yields an r ² of 0.49 (p-value = 0.04).

B_clim,a sensitivity to TP perturbations was not investigated by Sakai and Fujita (2017). Glaciers with a TP sensitivity < 1.0 m w.e. a⁻¹ ${20\,\%}^{-1}$ have either summer T ₂ < 0.5 °C (KYS, NAI, GIC, MZA, PIC, SIA) or summer TP ratios < 50 % (SOO, PL04, CSG, ABR). From the four glaciers ZHA, UG1, YAL and HAL with high summer T ₂ and summer TP ratios > 50 %, ZHA and YAL have a TP sensitivity > 1.7 m w.e. a⁻¹ ${20\,\%}^{-1}$ and ΔT ₂ range < 20 K, while UG1 and HAL have a TP sensitivity < 1.35 m w.e. a⁻¹ ${20\,\%}^{-1}$ and a ΔT ₂ range > 24 K. A r ² of 0.42 (p-value = 0.08) results from the multiple regression between the three explanatory variables summer T ₂, ΔT ₂ range and summer TP ratio and the TP sensitivity.

The ratio between the B_clim,a sensitivity to 1 K T ₂ change (ΔB _T2) and 20 % TP change (ΔB _TP) describes the relationship between temperature and precipitation sensitivity of the studied glaciers in comparison to each other (Table 5). ZHA, YAL, BIC, NAI and PIC have ratios above 1.4. In contrast, MZA and SIA are more sensitive to a 20 % TP change than to a 1 K T ₂ change with ratios below 0.7. All other glaciers have values between 0.7 and 1.3.

The ZHA sensitivity of −2.96 m w.e. a⁻¹ K⁻¹ in this study is more sensitive than the modelled sensitivity of −1.3 m w.e. a⁻¹ K⁻¹ by Zhu and others (2018). The difference is even more extreme for TP with a value of 1.96 m w.e. a⁻¹ ${20\,\%}^{-1}$ compared to 0.52 m w.e. a⁻¹ ${20\,\%}^{-1}$ of Zhu and others (2018). For PL04 the results are reversed with a modelled response of −1.28 m w.e. a⁻¹ K⁻¹ of Zhu and others (2018) and −0.55 m w.e. a⁻¹ K⁻¹ in this study and 0.29 m w.e. a⁻¹ ${20\,\%}^{-1}$ of Zhu and others (2018) and 0.68 m w.e. a⁻¹ ${20\,\%}^{-1}$ modelled here.

The comparison of the relative differences in sensitivity in this study confirms that glaciers with a high summer precipitation ratio are more sensitive (e.g. Sakai and others, 2015; Sakai and Fujita, 2017), and that summer precipitation ratio together with summer temperature and temperature range can be used to explain a large share of the climate sensitivity of glaciers in HMA.

5.3.3 Elevation range and area

The number of glaciers with different topographic features in this study is too low to statistically differentiate between climate background and topographic settings such as elevation range and area. Nevertheless, area does not seem to have a large effect on the sensitivity. For example, the magnitude of the B_clim,a sensitivity of PIC is similar to the sensitivities of the glaciers with an area smaller than 3.7 km², while the PIC is approximately 100 times bigger.

However, it seems the elevation range of the glaciers impacts glacier sensitivity. The smaller the elevation range, the stronger the relative effect of a shift in ELA (Sugden and John, 1976). MZA, SIA, NAI and CSG are the glaciers with the highest elevation range. They show the lowest overall SC_sea values. Furthermore, SC_sea of NAI (elevation range: 1686 m) and PL04 (elevation range: 1293 m) are lower than the SC_sea of the closely located glaciers HAL (elevation range: 400 m) and PL94 (elevation range: 577 m), respectively. However, the B_clim,a response of the NAI to uniform T ₂ perturbations (SC_uni) is higher than the B_clim,a response of the HAL. The reason might be that NAI has a large flat accumulation area, of which a distinct part may become an ablation area under uniform warming. BIC has the highest B_clim,a sensitivity of all glaciers to temperature and precipitation perturbations. Besides its climatic setting, this is likely due to its low elevation range (295 m).

When the elevation range is added to the summer TP ratio to explain the T ₂ sensitivity in a multiple regression, r ² increases from 0.29 to 0.38 (p-value = 0.05). The same applies if elevation range is added to summer T ₂ for the explanation of the precipitation sensitivities, with the r ² increasing from 0.26 to 0.38 (p-value = 0.04). The r ² between elevation range and temperature sensitivity is 0.3 (p-value = 0.03) and between elevation range and precipitation sensitivity is 0.37 (p-value = 0.01).

6. Conclusion

We simulated the annual glacier-wide climatic mass balance (MB-years 2001–2018) of 16 glaciers in HMA using a glacier energy- and mass-balance model. The model was calibrated to a geodetic glacier mass balance dataset using total precipitation scaling factors for each glacier, as total precipitation is the forcing variable with the highest uncertainty. We used a uniform system and the same approach for all glaciers to study interannual variability and climate sensitivity, and to discuss similarities and differences between glaciers covering the major mountain ranges in HMA. Despite the uncertainties in the climate forcing variables and the model parameterisations, the results provide estimates of the variability and sensitivity of glacier mass balance to climate fluctuations in uniform ways. The results indicate that it is crucial to separate liquid and solid precipitation due to their influence on the albedo parameterisation and refreezing within the snowpack. Sub-daily resolution of all processes and glacier melt in general are important. The sensitivity to temperature and precipitation cannot be analysed individually since both variables strongly affect ablation and accumulation in combined ways.

Glaciers at the margin of HMA generally have a higher interannual variability, while glaciers located to the southeast experience higher mass turnover. The results show that ELAs lower in elevation from South to North. All glaciers are most sensitive to changes in temperatures in summer compared to the rest of the year. In the case of precipitation, it is the change in spring and summer or only summer that is dominant. Precipitation sensitivity does not necessarily follow monthly accumulated snowfall. Only one glacier (Chhota Shigri Glacier) is most sensitive to precipitation perturbations in winter. The results show that it makes a distinct difference if accumulated snowfall or total precipitation is considered for the classification of the principal accumulation season of a glacier.

Glaciers with high mass turnover, high summer temperature, low annual temperature range and high summer to annual precipitation ratio are more sensitive to temperature changes, which confirms previous studies with simpler models on mass balance sensitivity in HMA. The sensitivity experiments in this study reveal non-linear relationships between climatic mass-balance response and both air temperature and precipitation perturbations. The ratio of refreezing to meltwater and rain fed into the snowpack falls for most of the glaciers between 29 % and 43 %. Glaciers with an elevation range > 1500 m and accumulation areas above 7000 m a.s.l. show higher ratios of > 51 %.

The resulting spatial patterns of mass-balance sensitivity and variability should be further investigated and considered in future glacier projections. Moreover, despite the high uncertainties and possible error compensation in this approach, the study demonstrates the possibility of such a uniform approach. The interactive coupling of COSIPY in HMA with higher resolution atmospheric datasets, including schemes to account for snowdrift and supraglacial debris cover could be the next steps to improve the approach. Within such an integrated setup, local measurements of any kind could be used in the future to motivate the usage of adapted values for the various model parameterisation. The approach offers the possibility to improve the representation of both large-scale and small-scale processes within a single and consistent physical framework. Furthermore, similar studies should be extended to hundreds of glaciers in the future to then statistically distinguish the differences between the influence of morpho-topographic features and climate forcing.