1. Introduction
The glaciers of High Mountain Asia (HMA) are important long-term hydrological reservoirs for municipal supply, irrigation, and hydropower (Pritchard, Reference Pritchard2019; Immerzeel and others, Reference Immerzeel2020; Li and others, Reference Li2022). Furthermore, they pose a significant threat through glacial lake outburst floods (GLOFs) and glacier involved rock and ice avalanches (Carrivick and Tweed, Reference Carrivick and Tweed2016; Furian and others, Reference Furian, Loibl and Schneider2021; Kropáček and others, Reference Kropáček, Vilímek and Mehrishi2021; Zheng and others, Reference Zheng2021).
Overall, HMA glaciers have been losing mass in recent decades (e.g. Neckel and others, Reference Neckel, Kropáček, Bolch and Hochschild2014; Kang and others, Reference Kang2015; Ke and others, Reference Ke, Ding, Li and Qiu2017; Shean and others, Reference Shean2020). Shean and others (Reference Shean2020) estimate a HMA-wide annual negative geodetic mass balance of (−0.19 ± 0.03) m w.e. a−1 between 2000 and 2018. Rapid warming is enhancing glacier retreat (Yao and others, Reference Yao2019). The warming rate of 0.34 K (decade)−1 between 1961 and 2010 (Wei and Fang, Reference Wei and Fang2013) 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, Reference Rounce, Hock and Shean2020). 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, Reference Miles2021).
Nevertheless, significant differences between the mass balances in different regions of HMA exist (e.g. Yao and others, Reference Yao2012; Li and others, Reference Li, Yao, Yang, Yu and Zhu2018; Bonekamp and others, Reference Bonekamp, de Kok, Collier and Immerzeel2019; Bhattacharya and others, Reference Bhattacharya2021). 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, Reference Hewitt2005; Farinotti and others, Reference Farinotti, Immerzeel, de Kok, Quincey and Dehecq2020), which is also visible in parts of the Pamir and West Kunlun Shan (Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; de Kok and others, Reference de Kok, Kraaijenbrink, Tuinenburg, Bonekamp and Immerzeel2020; Zhu and others, Reference Zhu2022). 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, Reference Yao2012; Mölg and others, Reference Mölg, Maussion and Scherer2014; Fugger and others, Reference Fugger2022), local circulation systems (Xu and others, Reference Xu2014; Curio and others, Reference Curio, Maussion and Scherer2015; Ma and others, Reference Ma, Lu, Chen, Pan and Hong2018), differences in snow to total precipitation ratios and precipitation seasonality (Bonekamp and others, Reference Bonekamp, de Kok, Collier and Immerzeel2019), and different glacier morphology like slope, aspect, topography and size (Yu and others, Reference Yu2013; Brun and others, Reference Brun2019).
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, Reference Yao2012). Generally speaking, the ISM influences the spring- and summer-accumulation type glaciers (Benn and Owen, Reference Benn and Owen1998; Fujita and Ageta, Reference Fujita and Ageta2000; Yao and others, Reference Yao2012). 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, Reference Yao2012).
Glaciers in wetter climates are more sensitive to temperature change (Oerlemans, Reference Oerlemans1997) and less sensitive to changes in precipitation (Yang and others, Reference Yang2013). Ablation and temperature are non-linearly connected (Bolibar and others, Reference Bolibar, Rabatel, Gouttevin, Zekollari and Galiez2022). 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, Reference Oerlemans and Knap1998; Fujita, Reference Fujita2008b). 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, Reference Qu and Hall2007; Fujita, Reference Fujita2008a; Ghatak and others, Reference Ghatak, Sinsky and Miller2014; Johnson and Rupper, Reference Johnson and Rupper2020; Arndt and others, Reference Arndt, Scherer and Schneider2021). 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, Reference Fujita2008b). Glaciers in arid, cold climates are less sensitive to temperature changes (Ohmura and others, Reference Ohmura, Kasser and Funk1992; Braithwaite, Reference Braithwaite2008). Glaciers with large winter precipitation are less sensitive to changes in the liquid/solid fraction of precipitation (Kapnick and others, Reference Kapnick, Delworth, Ashfaq, Malyshev and Milly2014). Such glaciers, as for example in the Pamirs, are more controlled by precipitation than by air temperature changes (Zhu and others, Reference Zhu2020).
Recent glacier mass balance sensitivity studies in HMA, confirm the higher sensitivity of summer accumulating glaciers (Sakai and others, Reference Sakai2015), a nonlinear temperature but linear precipitation sensitivity of glaciers (Wang and others, Reference Wang, Liu, Shangguan, Radić and Zhang2019), and that spatial differences in glacier response have a higher explanatory power than regional differences in climate change (Sakai and Fujita, Reference Sakai and Fujita2017). 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, Reference Sakai and Fujita2017).
These past studies, while informative, have been based on degree-day models (Wang and others, Reference Wang, Liu, Shangguan, Radić and Zhang2019), often at daily or coarser temporal resolutions, or focussing on the temperature sensitivity only (Sakai and Fujita, Reference Sakai and Fujita2017). Instead, we apply the physically-based COupled Snowpack and Ice surface energy and mass balance model in PYthon (COSIPY, Sauter and others, Reference Sauter, Arndt and Schneider2020) 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, Reference Muñoz-Sabater2021) 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 (Reference Shean2020) 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 (Reference Shean2020). 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 (Reference Maussion2014). We selected glaciers of different sizes because they are well-studied benchmark glaciers for specific regions. Shi and Liu (Reference Shi and Liu2000) 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, Reference Shean2020).
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 (Reference Shean2020) 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 (Reference Shean2020) for the period 2000 to 2018.
Full glacier names are provided in the caption of Fig. 1.
Table 2. ERA5-L mean surface pressure psfc, air temperature at 2 m T 2, relative humidity at 2 m RH 2, incoming shortwave radiation QSWin, incoming longwave radiation QLWin, wind speed at 2 m U 2, 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 (Reference Shean2020) 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 (Reference Shean2020).
The glaciers are Batysh Sook (SOO, Kenzhebaev and others, Reference Kenzhebaev2017), Zhadang (ZHA, Huintjes and others, Reference Huintjes2015), Urumqi Glacier No. 1 (UG1, Li and others, Reference Li2021), Yala (YAL, Stigter and others, Reference Stigter2018), Halji (HAL, Arndt and others, Reference Arndt, Scherer and Schneider2021), Bayi Ice Cap (BIC, Qing and others, Reference Qing, Han and Liu2018), Parlung No. 94 (PL94, Shaw and others, Reference Shaw2021), Keli Yanghe source (KYS), Naimona'nyi (NAI, Zhu and others, Reference Zhu2021), Parlung No. 4 (PL04, Shaw and others, Reference Shaw2021), Chhota Shigri (CSG, Azam and others, Reference Azam2016), Abramov (ABR, Denzinger and others, Reference Denzinger2021), Guliya Ice Cap (GIC, Thompson and others, Reference Thompson2018), Muztagh Ata (MZA, Holzer and others, Reference Holzer2015), Purogangri Ice Cap (PIC, Liu and others, Reference Liu2019) and Siachen (SIA, Berthier and Brun, Reference Berthier and Brun2019). 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 (Reference Shi and Liu2000). Eight glaciers have an area smaller than 4 km2, 12 glaciers smaller than 22 km2 and four glaciers are larger than 100 km2 (RGI6 RGI Consortium, 2017). According to the Shuttle Radar Topography Mission (SRTM, Farr and others, Reference Farr2007) 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−1 and KYS the most positive with +0.43 m w.e. a−1 (Shean and others, Reference Shean2020). 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, Reference Sauter, Arndt and Schneider2020). The surface energy balance (SEB) is solved with an optimisation algorithm and coupled through the surface temperature Ts, surface melt Msfc, 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:
where QM is the available melt energy, QSWin is the incoming shortwave radiation, α is the albedo, QLWin is the incoming longwave radiation, QLWout is the outgoing longwave radiation, QH is the sensible heat flux, QE is the latent heat flux, QG is the glacier heat flux and QR is the sensible heat flux of rain. All terms on the right-hand side of Eq. (1) are input (QSWin and QLWin) or intermediate (α) variables, or parameterised (QLWout, QH, QE, QG and QR) using the optimisation algorithm. The residual of the optimisation algorithm results in QM only if Ts 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 (Reference Cogley2011):
where csfc is surface accumulation, asfc is surface ablation, ci is refreezing, ai is subsurface melt, bsfc is the surface mass balance and bi is the internal mass balance. In COSIPY, csfc results from accumulated snowfall SFc 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 asfc. Refreezing of percolating water results in ci, and subsurface melt Msub resulting from penetrating radiation equals ai.
All types of frozen precipitation equal SFc 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 2 using a logistic transfer function (Hantel and others, Reference Hantel, Ehrendorfer and Haslinger2000). 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−2. 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 Bclim is the sum of all point bclim 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 Bclim,a and the glacier-wide cumulative climatic mass balance Bclim,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, Reference Cogley2011). 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 (Reference Oerlemans and Knap1998). Values for fresh snow, firn and ice albedo correspond to the value in Sauter and others (Reference Sauter, Arndt and Schneider2020). We apply the parameters from Mölg and others (Reference Mölg, Maussion, Yang and Scherer2012) for the parameters of time and height. For a comprehensive description of the model, please refer to Sauter and others (Reference Sauter, Arndt and Schneider2020).
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, Reference Muñoz-Sabater2021). From the forcing variables required for COSIPY (Table 3), only T 2 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 (Reference Arndt, Scherer and Schneider2021).
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 psfc to the mean elevation of the glacier as well as the interpolation to the GGPs. T 2 and T d,2 (intermediate variable to calculate relative humidity at 2 m RH 2) 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 2 to all GGPs, while for RH 2 no interpolation approach is applied. No downscaling is applied to incoming shortwave radiation QSWin and incoming longwave radiation QLWin. Also no interpolation is applied to QLWin. The radiation model after Wohlfahrt and others (Reference Wohlfahrt2016) is used. It corrects QSWin considering slope and aspect for the interpolation of QSWin to the GGPs. The scaling of TP is used as calibration of the Bclim,a to the remote sensing derived geodetic glacier mass balance of Shean and others (Reference Shean2020). The derived wind speed at 2 m U 2 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, Reference Thiel2020; Arndt and others, Reference Arndt, Scherer and Schneider2021) 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, SFc differs at each GGP, since TP is partitioned into wet and solid precipitation depending on T 2, 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, Reference Farr2007) 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, Reference Farr2007) 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, Reference Jarvis, Reuter, Nelson and Guevara2008) 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−1 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, Reference Arndt, Scherer and Schneider2021). 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−3 and a bottom snow density of 500 kg m−3. 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 2 and U 2 in accordance with Vionnet and others (Reference Vionnet2012). For the initialisation of the density and temperature profiles, please refer to the values given by Sauter and others (Reference Sauter, Arndt and Schneider2020). The resulting Bclim,a are compared with the geodetic mass balances of Shean and others (Reference Shean2020). If a glacier consists of multiple RGI6 outlines, results of Shean and others (Reference Shean2020) 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−1. 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 SCsea following Oerlemans and Reichert (Reference Oerlemans and Reichert2000) and uniform sensitivity characteristic SCuni similar to Schuler and others (Reference Schuler2005). For the SCsea, we first applied T 2 perturbations for each glacier to obtain an equilibrium Bclim,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, Reference Oerlemans and Reichert2000) 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 Bclim,a to the monthly perturbations.
For the SCuni, 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 SCsea, the reference runs for the SCuni are the simulations calibrated to Shean and others (Reference Shean2020).
4. Results
4.1 Simulated energy and mass balance
The simulated annual glacier-wide climatic mass balance Bclim,a, glacier-wide cumulative climatic mass balance Bclim,cum, annual accumulated snowfall SFc, 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 Bclim,a, PIC, PL04, PL94, ZHA, GIC have only negative Bclim,a and the rest have positive and negative Bclim,a. The glaciers with the most negative Bclim,cum (except ZHA) are glaciers that have positive and negative Bclim,a, in contrast to GIC and PIC on the Tibetan Plateau which have only negative Bclim,a but have not the most negative Bclim,cum over the whole study period. Bclim,a of GIC, PIC, ZHA and KYS – all located centrally within HMA – have the smallest interannual variability according to the coefficient of variation cv, with values < 0.15 (see Fig. 2). The detrended standard deviation of Bclim,a depends on the amount of Bclim,a. Therefore, cv is defined as the ratio between the detrended standard deviation of Bclim,a and the mean of the annual values of the mass turnover similar to Maussion and others (Reference Maussion2014). Bclim,a of glaciers in the Tian Shan occur with a high interannual variability with cv > 0.40. Bclim,a of the glaciers in the Himalayas are also highly variable with a mean cv > 0.35 including the highest interannual variability of NAI. Bclim,a of SIA also has a high interannual variability (cv = 0.44), while the Bclim,a glaciers in the Pamir, Pamir-Alay and Qilian Shan have a medium interannual variability (cv = 0.28, 0.29, 0.31). Bclim,a of Parlung glaciers have quite different interannual variability with cv values of 0.25 and 0.39.
Figure 2. Annual glacier-wide climatic mass balance Bclim,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 Bclim,a and coefficient of variation cv (σ of Bclim,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 Bclim,cum, while the ZHA with the most negative Bclim,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 SFc, 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−1 and a mean value (MB-year 2001–2018) of only 0.007 m w.e. month−1 in July. Using SFc 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 QSWnet is the greatest energy source for all glaciers. GIC, and ZHA located at the interior of HMA receive the highest QSWnet compared to the other glaciers. The glaciers with the most positive Bclim,a KYS receive the lowest energy input from QSWnet because of high glacier-wide albedo. Furthermore, glaciers in the same mountain ranges receive energy from QSWnet in the same range. Glaciers in the Himalayas receive between 57 W m−2 and 64 W m−2, the Parlung glaciers between 47 W m−2 and 54 W m−2, and glaciers in the Tian Shan, Pamir and Pamir-Alay between 44 W m−2 and 51 W m−2. Net longwave radiation QLWnet 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 QE of each glacier is an energy sink, while QH is an energy source for some glaciers and an energy sink for the others. The sign of QE and QH 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, QG is always positive and directed to the surface in the overall mean. QR plays a minor role and is certainly always positive. QSWnet, QLWnet and the closing term of the energy balance QM 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 Bclim,a components. The most important mass flux is Msfc for all glaciers with the exemption of KYS, which has the least negative annual rate with −0.44 m w.e. a−1. The most negative mean annual value is YAL with −2.83 m w.e. a−1. With the exemption of MZA, refreezing is higher than SFc, while it has to be kept in mind that refreezing depends on the available percolation water through Msfc, 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 SFc. Msub and sublimation have a minor influence on the mass budget. Despite high QE 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 106 J kg −1 K −1 $\gg$
latent heat of fusion: 3.337 105 J kg −1 K −1).
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 bclim,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−1, while the bands with the highest ablation rates reach 6 m w.e. a−1. 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 Bclim,a, SFc, 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 Bclim,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 SCsea 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 Bclim,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 (Reference Oerlemans and Reichert2000). Red bars are the dependence of the annual glacier-wide climatic mass balance Bclim,a on monthly temperature perturbations of 1 K and blue the dependence of Bclim,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 SCsea 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 SFc (see Fig. 3). Some differences are evident, though, such as SFc 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 SFc is very low in these months. PIC, BIC and UG1 show a peak of SFc and of precipitation SCsea 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 SFc is in February.
Similar to Schuler and others (Reference Schuler2005) the SCuni are shown in Fig. 7. Noteworthy is the shape of the SCuni 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 2, 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 2 perturbations. Noteworthy is the Bclim,a response of nearly −8 m w.e. a−1 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 SCsea to reach a zero Bclim,cum (see Fig. S7) and not the ones scaled to Shean and others (Reference Shean2020).
Figure 7. Uniform sensitivity characteristics SCuni. Displayed are the annual glacier-wide climatic mass balance Bclim,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 2 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, Reference Shaw2021). The error in QSWin 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 QLWin output is used directly, but better modelling of cloud cover and QLWin is an important objective (Kok and others, Reference Kok2019). Scaling of U 2 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, Reference Immerzeel, Wanders, Lutz, Shea and Bierkens2015; Gao and others, Reference Gao, Xiao, Chen, Xu and Zhang2018; Hamm and others, Reference Hamm2020; Arndt and others, Reference Arndt, Scherer and Schneider2021). To account for orographic effects on precipitation, high resolution precipitation modelling (Hamm and others, Reference Hamm2020) or at least information about region-specific lapse rates (Machguth and others, Reference Machguth, Paul, Kotlarski and Hoelzle2009) would be needed. A larger amount of SFc is simulated in the higher-elevated areas compared to the lower-elevated areas of the glaciers due to the lapse rate of T 2 and the fact that T 2 separates TP into either rain or SFc.
Of the COSIPY subroutines, Bclim,a is most sensitive to the albedo parameterisation (e.g. Arndt and others, Reference Arndt, Scherer and Schneider2021), as is typical for such models (e.g. Zhu and others, Reference Zhu2018). Different climate regimes and SFc and T 2 patterns lead to different constants for the albedo parameterisation (Oerlemans and Knap, Reference Oerlemans and Knap1998) for individual glaciers. Nevertheless, due to limited data it is not possible to use glacier-specific parameters. Bclim,a is sensitive to the separation of TP into solid and liquid based on the logistic transfer function, which is highly dependent on T 2 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, Reference Bintanja and Reijmer2001; Zhu and others, Reference Zhu2018). Debris-cover enhances melting (lower albedo) until a certain thickness and reduces melting (isolating effect) above the threshold (Östrem, Reference Östrem1959). Further, it changes the energy balance because surface temperatures above the melting point temperature are possible (Fugger and others, Reference Fugger2022). SIA is covered between 3 % and 13 % with debris (Agarwal and others, Reference Agarwal2017; Kumar and others, Reference Kumar, Negi and Kumar2020). 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, Reference Laha2017) 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, Reference Houska, Kraft, Chamorro-Chavez and Breuer2015) was used to quantify the global sensitivity of Bclim,a of the 16 glaciers based on the Fourier amplitude sensitivity test (FAST) of Saltelli and others (Reference Saltelli, Tarantola and Chan1999). 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 Bclim,a deviation of −0.11 m w.e. a−1 (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−1 and 1.94 m w.e. a−1 for climate forcing variables and between 0.46 m w.e. a−1 and 2.29 m w.e. a−1 for model parameters. Zhu and others (Reference Zhu2018) calculated a change in mass balance of 0.92 m w.e. a−1 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 (Reference Zhu2018) regarding the extremely high sensitivity of ZHA compared to PL94. Wang and others (Reference Wang, Liu, Shangguan, Radić and Zhang2019) calculated a mass balance response (averaged over 45 glaciers) of −1.59 m w.e. a−1 from varying temperature lapse rates by −20 % and the least response of −0.45 m w.e. a−1 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 (Reference Wang, Liu, Shangguan, Radić and Zhang2019) with monthly resolution and in the case of Zhu and others (Reference Zhu2018) 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 Bclim,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−1. 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−1. 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 Bclim,a could reduce the resulting equifinality problem (Barandun and others, Reference Barandun2018).
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 2 because U 2 is underestimated in the dataset. A comparison between the mean, scaled U 2 and the mean measured U 2 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 (Reference Muhammad and Tian2020) measured a TP amount of 185 mm a−1 for one year compared to the 2000–2018 mean of 246 mm a−1 here. The scaled TP amount of YAL is 2079 mm a−1, which is higher than reported by nearby in situ observations (cf. Stumm and others, Reference Stumm, Joshi, Gurung and Silwal2021). However, taking into account the possible undercatch of snow (Rasmussen and others, Reference Rasmussen2012), the average value of 1867 mm a−1 (Sharma and others, Reference Sharma, Hamal, Khadka and Joshi2020, 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 (Reference Maussion2014).
The mean energy fluxes QSWnet, QLWnet, QH and QE for ZHA are all between the maximum and minimum values given in four other studies summarised in Table 7 in Zhu and others (Reference Zhu2018). Mölg and others (Reference Mölg, Maussion, Yang and Scherer2012) and Huintjes (Reference Huintjes2014) include shortwave penetrating radiation in the calculation of QG, in contrast to this study. Adjusting QG to include shortwave penetrating radiation in this study results in values for QG being similar to both studies and with the same sign. The resulting QM value of 21.9 W m−2 for ZHA is close to the mean value (19.4 W m−2) of the four other studies (maximum: 30.6 W m−2, minimum: 13.7 W m−2).
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 2 = 0.59 and r 2 = 0.4) by the modelling chain. For PL94, the result is not significant with r 2 = 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 Bclim,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, Reference Rounce, Hock and Shean2020) with multiple indicators (e.g. Barandun and others, Reference Barandun2018) was applied.
Figure S13a shows the comparison of the mean Bclim,a simulated over the entire study period with the estimated elevation changes of Hugonnet and others (Reference Hugonnet2021). It shows a good (r 2 = 0.68) agreement of the change in the different mountain ranges. However, the dataset of Hugonnet and others (Reference Hugonnet2021) is not strictly independent, since both the estimate of Shean and others (Reference Shean2020) and Hugonnet and others (Reference Hugonnet2021) are derived from ASTER DEMs. Furthermore, we compared simulated ELAs and ELAs derived by Miles and others (Reference Miles2021) of the studied glaciers in Fig S13b. Again, the significant correlation (r 2 = 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 (Reference Fujita and Ageta2000) 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 (Reference Li, Yao, Yang, Yu and Zhu2018) simulated a rate of 25 % for Qiangtang No.1 Glacier, which is also most similar to the PIC. Yang and others (Reference Yang2013) 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 (Reference Veldhuijsen2021) 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 (Reference Fujita and Ageta2000), Li and others (Reference Li, Yao, Yang, Yu and Zhu2018) and Yang and others (Reference Yang2013).
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 Bclim,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 (Reference Shean2020, 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, Reference Greuell and Smeets2001; Oerlemans and Klok, Reference Oerlemans and Klok2002). Further, the expectation of higher sublimation rates for more continental type glaciers (e.g. Zhu and others, Reference Zhu2015) 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 QG. A northwest–southeast pattern is most evident in Bclim,a and total mass turnover (see Fig. S6).
The definition of the accumulation season of the glaciers depends on whether TP or SFc is taken as a classifying variable. Some glaciers are classified differently here than in other studies since we refer to SFc only. Considering TP, most of the glaciers (except SIA, CSG, PL04, PL94) are summer-accumulation type glaciers, but the pattern changes considering only SFc. For example, HAL, ZHA and PIC are not clear summer-accumulation type glaciers using SFc 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 (Reference Yang2013) and Jouberton and others (Reference Jouberton2022). The percentage of TP occurring as SFc in this study fits the pattern of Maussion and others (Reference Maussion2014), 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 (Reference Oerlemans and Reichert2000). The results point towards higher overall T 2 sensitivity of ABR compared to Oerlemans and Reichert (Reference Oerlemans and Reichert2000). 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 (Reference Oerlemans and Reichert2000), T 2 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 2 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 (Reference Oerlemans and Reichert2000). Further, sensitivity increases when the climate is wetter (Oerlemans, Reference Oerlemans1997; Fujita, Reference Fujita2008a) as a general pattern. Moreover, SCsea reveals that Bclim,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 SFc in summer is present. All glaciers with a high summer precipitation amount are very sensitive to summer changes in T 2, which is in accordance with Naito (Reference Naito2011). Total mass turnover is important for the SCsea as well, which is supported by the statistically significant correlation (r 2 = 0.38, p-value = 0.01) between the mass turnover and overall SCsea (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 Bclim,cum at KYS. Only in August and to some extent June, T 2 temperature perturbations have an influence on Bclim,a.
5.3.2 Uniform sensitivity characteristics
Figure 7 reveals that neither the sensitivity to T 2 nor to TP is linear (cf. Oerlemans and Fortuin, Reference Oerlemans and Fortuin1992; Fujita, Reference Fujita2008a). This partially contradicts Wang and others (Reference Wang, Liu, Shangguan, Radić and Zhang2019), 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 (Reference Sakai and Fujita2017) identified summer temperature, temperature range and summer precipitation ratio as the variables with the most explanatory power for the Bclim,a sensitivity to T 2. The Bclim,a response of glaciers with a summer TP ratio < 50 % (SOO, PL94, PL04, CSG, ABR, MZA, SIA) results in a low T 2 response. ZHA, UG1, YAL and HAL have summer precipitation ratios > 50 %. Nevertheless, from the four glaciers UG1 and HAL have high ΔT 2 ranges resulting in a low T 2 response, while ZHA and YAL have low ΔT 2 ranges resulting in a high T 2 response < −2.95 m w.e. a−1 K−1.
Table 5. Topographic and climatic characteristics and sensitivity indices.
Elevation range (ER); summer temperature (T 2 JJA) calculated from the long-term (MB-year 2001–2018) mean of June, July and August; temperature range (ΔT 2) 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 2 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 SCuni. 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 2 of −5 °C and −0.9 °C, respectively, resulting in low and medium T 2 sensitivity. BIC has the highest T 2 response with a value of −3.4 m w.e. a−1 K−1. Only high summer T 2 and summer precipitation ratio > 50 % in combination can explain this high sensitivity to some extent. KYS has a low summer T 2 resulting in the very low sensitivity even though the summer TP ratio is 60 %. The influence of summer TP ratio on T 2 sensitivity is supported by an r 2 of 0.29 (p-value = 0.03). A multiple regression model between the three explanatory variables summer T 2, ΔT 2 range and summer TP ratio and the T 2 sensitivity yields an r 2 of 0.49 (p-value = 0.04).
Bclim,a sensitivity to TP perturbations was not investigated by Sakai and Fujita (Reference Sakai and Fujita2017). Glaciers with a TP sensitivity < 1.0 m w.e. a−1 ${20\,\%}^{-1}$
have either summer T 2 < 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 2 and summer TP ratios > 50 %, ZHA and YAL have a TP sensitivity > 1.7 m w.e. a−1 ${20\,\%}^{-1}$
and ΔT 2 range < 20 K, while UG1 and HAL have a TP sensitivity < 1.35 m w.e. a−1 ${20\,\%}^{-1}$
and a ΔT 2 range > 24 K. A r 2 of 0.42 (p-value = 0.08) results from the multiple regression between the three explanatory variables summer T 2, ΔT 2 range and summer TP ratio and the TP sensitivity.
The ratio between the Bclim,a sensitivity to 1 K T 2 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 2 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−1 K−1 in this study is more sensitive than the modelled sensitivity of −1.3 m w.e. a−1 K−1 by Zhu and others (Reference Zhu2018). The difference is even more extreme for TP with a value of 1.96 m w.e. a−1 ${20\,\%}^{-1}$
compared to 0.52 m w.e. a−1 ${20\,\%}^{-1}$
of Zhu and others (Reference Zhu2018). For PL04 the results are reversed with a modelled response of −1.28 m w.e. a−1 K−1 of Zhu and others (Reference Zhu2018) and −0.55 m w.e. a−1 K−1 in this study and 0.29 m w.e. a−1 ${20\,\%}^{-1}$
of Zhu and others (Reference Zhu2018) and 0.68 m w.e. a−1 ${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, Reference Sakai2015; Sakai and Fujita, Reference Sakai and Fujita2017), 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 Bclim,a sensitivity of PIC is similar to the sensitivities of the glaciers with an area smaller than 3.7 km2, 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, Reference Sugden and John1976). MZA, SIA, NAI and CSG are the glaciers with the highest elevation range. They show the lowest overall SCsea values. Furthermore, SCsea of NAI (elevation range: 1686 m) and PL04 (elevation range: 1293 m) are lower than the SCsea of the closely located glaciers HAL (elevation range: 400 m) and PL94 (elevation range: 577 m), respectively. However, the Bclim,a response of the NAI to uniform T 2 perturbations (SCuni) is higher than the Bclim,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 Bclim,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 2 sensitivity in a multiple regression, r 2 increases from 0.29 to 0.38 (p-value = 0.05). The same applies if elevation range is added to summer T 2 for the explanation of the precipitation sensitivities, with the r 2 increasing from 0.26 to 0.38 (p-value = 0.04). The r 2 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2023.46
Data
The data presented in this study are available on request from the corresponding author.
Acknowledgements
This research was financed by the German Research Foundation's (DFG) research grants ‘Precipitation patterns, snow and glacier response in High Asia and their variability on sub-decadal time scales, sub-project: snow cover and glacier energy and mass balance variability’ (prime-SG, SCHN 680/13-1), ‘Dynamic Response of Glaciers in the Qilian Shan to Climate Change’ (Dyn-Q, SCHN 680/17-1), and ‘Glacial lake outburst floods in the Halji region, Nepal’ (Halji, SCHN 680/19-1). We use data published by Shean and others (Reference Shean2020), who are thanked for sharing their data online. We are grateful to Evan Miles, Scientific Editor, Nicolas Cullen, Associate Chief Editor and three anonymous reviewers for their extensive, insightful and very constructive feedback that helped to substantially improve the manuscript in terms of scientific content, structure and language.
Authors’ contributions
Both authors jointly developed the concept and methodology of this study. The same accounts for the discussion and interpretation of the results. AA collected and prepared the data, executed the simulations and created the plots and tables. AA wrote the first draft of the manuscript. CS edited and reviewed the manuscript. CS acquired the funding and supervised the projects. Both authors have read and agreed to the published version of the manuscript.
Open access
