Introduction

In High Mountain Asia (HMA), present and future glacier down wasting and retreat will lead to changes in seasonal and spatial patterns of meltwater contributions to streamflow1,2. In highly glacierised catchments, changes in meltwater inputs will considerably affect the timing and magnitude of river discharge. Given the large population densities in the region and the strong socio-economic dependence on streamflow for irrigation, industry, and power generation, reliable estimates of future glacier mass change are critical.

Ablation represents the ensemble of processes that lead  to ice and snow mass loss. It includes melt, sublimation or evaporation at the surface, as well as wind-driven transport and sublimation of blowing snow3. By implicitly assuming that surface melt is the dominant ablation term, hydrological and glaciological models often use temperature-index (TI) or enhanced temperature-index (ETI) approaches to quantify total ablation. Such models conveniently rely only on near-surface air temperature, as well as a degree-day factor (TFTI) which depends on the state of the surface i.e., ice or snow, and assume that melt is a linear function of the sum of the difference between the daily mean air temperature Ta and an air temperature threshold (TTI)4,5. This relationship is generally well-established at lower altitudes. The ETI model builds on the TI approach by including the contribution of the daily mean net incoming shortwave radiation (SWnet) multiplied by a shortwave radiation factor (SRF, kg day−1 W−1 6). At high elevation glaciers, such as those found throughout High Mountain Asia (HMA), incoming shortwave radiation (SWinc) is the dominant energy input7 and the SEB relates only partly to daily mean Ta, so numerous studies use ETI in place of TI. Estimates of albedo (α) and SWinc are required for an ETI approach, and  SWinc can be approximated  as theoretical potential shortwave incoming at the top of the atmosphere. Atmospheric transmissivity, cloudiness and α can be estimated from empirical models calibrated on Ta and relative humidity (RH) data (see Methods for a detailed description).

However, surface melt is more suitably modeled with the net surface energy balance (SEB), which accounts for the energy entering the ice or snow surface from above (downward, positive) and that leaving it (negative), upwards to the air, or downwards to the deeper layers, by heat conduction and transport. TI and ETI calibrations depend on the driving components of the SEB. These differ over time, latitude, and elevation. At mid-latitude and low-elevation sites, the TI is reliable since SWinc, net longwave radiation (LWnet) and the turbulent flux of sensible heat (H) are well correlated with Ta and have a strong control on glacial and snow melt. At high elevation sites, the links between glacial melt and Ta or SW and the TI (or ETI) relations have rarely been validated8. Studies of glacier SEB in the HMA7,9,10,11,12,13,14, and other high-altitude regions15,16,17 remain scarce (Table 1). On top of net shortwave radiation (SWnet), they demonstrate the importance of LWnet as a loss and indicate relatively small turbulent heat fluxes, which are generally positive for sensible heat (H) and negative for latent heat (LE).

Table 1 Daily mean contributions from net shortwave and longwave radiation and turbulent fluxes from various energy balance studies in High Mountain Asia.
Full size table

Sublimation, wind-driven snow redistribution, and blowing snow sublimation processes must be accounted for as part of net ablation3,18. Such processes have previously been shown to be important for glacier ablation in the region18,19. These processes are not directly linked to Ta or SW, and calibrations of TI and ETI upon observed ablation are not necessarily consistent in time and space when melt is not the only important contributor4,20,21. The numerical simulation of these additional ablation processes is necessary and possible through physical models but requires knowledge of many physical variables. Their use is thus severely constrained by in-situ data availability. That is why the TI and ETI models are preferably employed in HMA22,23 to quantify melt or even total ablation. Here, we question their validity, and we provide guidelines for future use of these models in the region. Our study relies on in-situ measurements of glacier surface lowering and energy balance obtained in the Nepalese Himalaya (Fig. 1). In the Results, we describe our observations of surface height and α changes, surface temperature (Ts) on the glacier, and the other surface characteristics (Fig. 2), together with the observed meteorological conditions (Fig. 3). We then show statistics for the SEB and its components (Fig. 4). Subsequently we present correlations between daily SEB, its components and on-glacier Ta or α (Fig. 5, see also Table S2). This permits us to assess the significance of key individual meteorological and surface variables in the SEB, and the relevance of TI and ETI parameters for surface ablation at different sites and for the different seasons. In the Discussion, we calibrate TI and ETI models to the observed surface ablation as derived from glacier surface height measurements and discuss the temporal and spatial consistency of their parameters (Fig. 6, see Methods for details). We analyse the statistical relations between Ta, SEB, glacial melt and ablation in terms of physical processes. This provides a critical view of ablation modelling with TI and ETI models in HMA environments. In the Conclusions we propose guidelines for future studies. Details about the Methods are provided in the final section.

Figure 1
figure 1

Study area map with elevation (Z), location of measurement sites (circles), and RGI5.0 glacier extents (light blue).

Full size image
Figure 2
figure 2

Summary of meteorological and surface lowering data from the two ablation zone sites (Yala Glacier, 5350 m a.s.l. and Mera Glacier, 5380 m a.s.l.) and the one accumulation zone site (Mera Glacier, 6352 m a.s.l.) between 2013 and 2017.

Full size image
Figure 3
figure 3

Daily averages of hourly measurements from the meteorological variables at the three on-glacier automatic weather stations: Mean incoming shortwave (black) and longwave (orange) radiation, mean air temperature (red), mean wind speed (purple) and relative humidity (blue). Mera Glacier, 6352 m a.s.l., was buried in snow on 28th of July 2016 resulting into extreme or no average variables.

Full size image
Figure 4
figure 4

Boxplot summaries of mean daily surface energy balance components calculated from all available data, classified by season and grouped into ablation (upper panel) and accumulation (lower panel) sites. Each box upper (resp. lower) edge is drawn at the third quartile (resp. first) of the variable distribution, and whiskers provide the same value plus (resp. minus) the interquartile range. Values outside the whiskers are provided by the dots.

Full size image
Figure 5
figure 5

Values of the Pearson’s correlation coefficient between the daily mean SEB and daily mean air temperature (R(Ta, SEB)), the daily mean SEB and daily mean SEB components (R(x, SEB)) and the daily mean SEB and α (R(α, SEB)), during the three defined seasons. Each bar represents a season’s value for one given year. The circles joined with the black line indicate the value of the correlation coefficients obtained when combining the data from all the years for a given season and site.

Full size image
Figure 6
figure 6

Cumulated ablation calculated with the surface lowering measurements (thick blue line), with the surface energy balance for changing z0 values (orange dashed and continuous lines), with the TI (red line) and ETI (clear blue line) with one fixed set of factors. The hourly wind speed is shown upside down (green curve). Periods of surface melt (Ts = 0) are highlighted in orange. Results from Mera Glacier, 5380 m a.s.l in 2014 and 2017 (a) from Yala Glacier, 5350 m a.s.l., in 2014, 2016 (b) and Mera Glacier, 6352 m a.s.l, in 2015 and 2016 (c).

Full size image

Results

Observed climatology and surface height changes

The Nepalese Himalayas are heavily influenced by the South Asian monsoon24, which brings most of the annual precipitation to the region between June and September. Following this annual pattern we defined 3 seasons with identical start and end dates each year. This is optimized for data availability, maximizing duration of overlapping time ranges with simultaneous observations of surface characteristics and meteorology (Figs 1 and 2). The pre-monsoon period, characterized by increasing humidity and temperatures23 starts 10 May and finishes 5 June. The monsoon is defined as 6 June–30 September25, and the post-monsoon period extends from 1 October to 15 November. We do not discuss the winter season in detail. Data availability and measurement details are found in Table 2.

Table 2 Summary of meteorological measurements used in this study.
Full size table

Annual surface lowering of up to 2 m was observed at the lowest ablation zone sites of Yala and Mera Glaciers (Fig. 2). In the ablation zone, this height change can be attributed to mass loss through mainly surface melt, as daily maximum surface temperatures (Ts) reach 0 °C (Fig. 2., see Methods for the calculation of Ts). However, surface lowering also occurs when daily maximum Ts is below 0 °C. This indicates that processes such as snow compaction, wind-driven erosion or sublimation must be occurring, as suggested by previous studies over glaciers in the HMA18,19,26. Also, surface lowering rates were variable in time (Fig. 2). Precipitation events resulted in fresh snow deposition, as indicated by increases in surface height: these events also increased α, and were more frequent during monsoon than during other seasons. During the monsoon, the high-elevation accumulation zone site on Mera Glacier experienced both surface height gains and losses. Surface melt is a possible contributor to this surface lowering, since the Ts occasionally reached 0 °C (Fig. 2). However, observations of densified ice layers inside the snowpack at this site indicate that water refroze in the lower layers: the glacier remains cold at this elevation18 and ablation therefore likely occurred through other processes.

The large scale atmospheric forcing was similar over both glaciers, as reflected by the meteorological measurements which evolved synchronously23 (Fig. 3). At all sites, incoming shortwave radiation (SWinc) was greatest at the beginning of June27, with an average daily value of 304 W m−2 for all pre-monsoon periods between 2010–201528. However, SWinc quickly decreased to an annual minimum at the summer solstice. This is due to increasing cloud cover starting from the end of May, as humid air masses are transported northwards by the South Asian monsoon, which is active until September in central Nepal24. In October, SWinc increased rapidly as the monsoon clouds and humidity retreated and clear skies returned to the region. SWinc then followed a standard seasonal pattern with a second minimum at the winter solstice. At the Mera accumulation zone site, SWinc was generally more intense than at the ablation zone sites, and this is likely due to Sun’s radiation   reaching ground after having crossed less atmospheric layers.

At all sites, high elevation and cold clear skies favor low values of incoming longwave radiation (LWinc), especially during winter17 (Fig. 3). Since clouds emit large amounts of thermal radiation29,30 cloudiness is strongly linked to LWinc at high altitudes. Cloudiness increased progressively through spring to reach a plateau in summer, as shown by the LWinc increase. At all sites, daily mean Ta rose synchronously with net longwave radiation (LWnet), starting from a minimum at the end of January and reaching a plateau of its highest values at the beginning of June each year. In the ablation zones of both glaciers daily mean Ta varied between −10 °C and 0 °C, but Ta rarely exceeded 0 °C as any excess energy was converted to snow or ice melt. Following the end of monsoon, Ta and LWinc decreased rapidly, though Ta showed large variations during autumn and winter. At Mera Glacier accumulation zone AWS (6352 m a.s.l), located about a thousand meters higher than the ablation zone sites (Fig. 1, Table 2), mean daily Ta values are consistently lower, and temperatures below −25 °C were observed in winter. The mean daily wind speed (w) remained below 10 m s−1 at 5350 m a.s.l. on Yala Glacier, while at 5380 m a.s.l on Mera Glacier, mean daily w greater than 15 m s−1 were observed. The accumulation zone site at Mera Glacier (6352 m a.s.l.) experienced the greatest mean daily w (>20 m s−1). At all sites the maximum wind speed was observed during winter and pre-monsoon, and is related to westerly disturbances and south to north shifts in the position of the jet stream. Daily mean w was the lowest during monsoon (~2 m s−1). Daily mean relative humidity (RH) was similar at all sites. Its variability was greatest during the pre-monsoon, and lowest during the monsoon. In all years, RH increased progressively from its minimum (~20%) in January to a plateau of 90–100% between June and September. RH then decreased abruptly in October, back to values around its minima.

Ablation Zone Surface Energy Balance and Mass Loss

Below, we present the statistics of the SEB components. The temporal correlations of total surface energy balance (SEB) with its components x and with Ta, α, are noted as R(x, SEB). Partial correlations are noted R(y, SEB)x (where y can be any variable, see Methods). The partial correlation of y with SEB indicates whether or not the correlation of y with SEB is due to a strong correlation of y with x. SEBs were strongly site- and season-dependent (Fig. 4). On average, pre-monsoon SEBs were positive at ablation zone sites. There the main energy input during pre-monsoon was +SWnet (daily mean ~+100 Wm−2). LWnet was a loss of energy of about −50 Wm−2, similar to that observed on tropical glaciers in South America17,31, in a similar climatic setting. Sensible turbulent fluxes (H) were positive, but latent heat fluxes (LE) were negative. On average, there was a slight net energy loss due to turbulent energy fluxes (H + LE ~ −15 Wm−2). Positive net SEBs and surface temperatures of 0 °C (Fig. 2) in the ablation zone suggest that the observed surface lowering was mainly driven by melt.

Variations in the SEB between sites and years were negatively correlated with LWnet (R(LWnet, SEB) = −0.52, Fig. 5, Table 3) and poorly or not at all correlated with the turbulent fluxes (R(H, SEB) = 0.12, R(LE, SEB) < 0.0). We assume that melt controls ablation when Ts = 0 °C, and directly convert positive SEB into surface melt. We obtain an average melt rate of −41 kg m−2 day−1 (Table 4). However, in the ablation zone mass loss also occurred through evaporation/sublimation, as indicated by the negative values of LE which correspond to an average mass loss of −3 kg m−2 day−1 over all the pre-monsoon periods (Table 4). During monsoon, the SEB increased (Figs 3, 4), with net daily energy gains greater than 100 Wm−2. Substantial increases in LWinc due to monsoon cloud cover were largely responsible for the SEB gains, as monsoon LWnet was only slightly negative. Surface ablation was large (−75 kg m−2 day−1 on average) and mainly due to melt as Ts was 0 °C. Turbulent fluxes were low or null. Differences in the SEB between the two ablation sites were likely due to varying SWnet (R(SWnet, SEB) = 0.90). This suggests that a fast alternation between snow and ice surfaces exerts a strong control on the SEB through α changes (R(α, SEB) = −0.80). Monsoon SEBs are weakly correlated with LWnet (R(LWnet, SEB) = 0.15), which signifies that cloudiness and higher longwave radiation totals are relatively stable. Nevertheless, the SWnet controls on the SEB are rather linked to changes in LWnet than to changes in α (R(SWnet, SEB)LWnet = 0.19 against R(SWnet, SEB)a < 0.01, Table S2). During post-monsoon at the ablation zone sites, clear sky conditions and high α (>0.8) after monsoon snowfalls reduce SWnet (often below 100 Wm−2), and large losses of energy occurred through large negative LWnet (<−50 Wm−2). The controls of LWinc over SEB are large32, reaching R(LWinc, SEB) = 0.66 at Yala ablation site (Table 3). Net turbulent fluxes were also negative, and dominated by negative latent heat fluxes, favored by intense wind. Small ablation totals were likely due to sublimation as Ts was mostly negative during the post-monsoon (Fig. 2), and net SEB was often negative.

Table 3 Correlation coefficients between surface energy balance components, air temperature and albedo, cumulating all available data.
Full size table
Table 4 Summary of the governing ablation processes and their relative contribution to the ablation (assuming no wind erosion).
Full size table

Accumulation Zone Surface Energy Balance and Mass Loss

The accumulation zone SEB is notably different from ablation zone sites during the main melt periods (pre-monsoon and monsoon). At the Mera Glacier accumulation zone site, SWnet was the main source of energy in the pre-monsoon. However, SEB remained negative throughout the pre-monsoon and surface melt was zero as Ts < 0 °C (Fig. 2). As this site (6352 m a.s.l.) is higher than the other sites, LWinc was lower and LWnet was more negative (LWnet ~ −100 Wm−2). Turbulent fluxes were more intense than at the lower elevation sites due to higher wind speeds, and RH values were generally lower (Fig. 3). As a result, net turbulent fluxes were negative (−25 Wm−2 on average), and reached minimum values of −80 Wm−2. Observed surface lowering was due mainly to wind-driven erosion, compaction, or sublimation. During the monsoon, SWnet at the Mera Glacier accumulation zone site was lower than during pre-monsoon. Reduced atmospheric transmissivity due to clouds and humidity resulted in decreased SWinc (50 to 100 Wm−2), while frequent fresh snow falls (Fig. 2) increased α and thus SWout. The cloud emission raised LWinc, but not as much as at the lower sites. H + LE was still negative and variable through the monsoon. Accumulation zone net SEB was thus lower than observed at ablation zone sites, and more variable. Our calculations indicate that sublimation as inferred from the latent heat fluxes (−12 kg m−2 day−1, Table 4) is a larger contributor to mass loss than melt (−9 kg m−2 day−1). While wind-driven ablation processes must be a significant factor for surface lowering and is not captured by our set up, enhanced sublimation rates for wind-blown snow18 suggest that our surface sublimation rate is a conservative estimate.

Air temperature, surface energy balance, melt and ablation

Correlations between Ta, SEB, and SEB components provide insights into the mechanisms driving the variation of the melt. At the ablation zone sites on Yala and Mera Glaciers during pre-monsoon, variations in the mean daily SEB were highly correlated to changes in mean daily SWnet (R(SWnet, SEB) ~0.6 to 0.9, Fig. 5), and less but still significantly, to daily mean Ta (R(Ta, SEB) ~0.3 to 0.5). The contribution of the statistical relation existing between Ta and SWnet to the overall R(Ta, SEB) (noted R(Ta, SEB)SWnet, Methods), was between 35% and 42% (Table S3). Mean daily Ta was highly correlated to SWnet as days with high insolation totals were warmer. This, in turn, enhances the entire SEB through increased LWinc and H. During the monsoon, Ta is similarly related to the SEB, but SWnet is better correlated to the SEB than during pre-monsoon. We observe a high negative correlation of α with the SEB at ablation zone sites (Fig. 5), for which daily mean Ta was a poor indicator. In the accumulation zone at Mera Glacier (6352 m a.s.l.), during pre-monsoon, Ta and SWnet variations related much less to the SEB than at the low elevation, ablation zone sites (Fig. 5). Overall, individual radiation terms did not explain the small SEB changes well. Similarly, during the monsoon, Ta as well as SWnet variations are less correlated to the SEB (R ~ 0.3) than at the lower elevation sites. Permanent snow cover at the high elevation site drives high α, and relatively low variability in SWnet. LWnet also exhibited a weak link with the SEB. However, the net radiation related well to the net SEB. Changes in the SEB related well to changes in H (R(H, SEB) ~0.49) and LE (R(LE, SEB) ~0.70), significantly more than in the ablation zone, as winds were more intense (Figs 3, 6). There was not enough data for a post-monsoon analysis.

Discussion

The performance of the ETI and TI using only one set of parameters for modelling the observed ablation at different sites and periods is limited (Fig. 6). The lack of consistency in TI and ETI parameters between sites and years (Table S4) is similarly problematic. For many cases, the calibration procedure even fails to optimize the parameters. Ablation modeled with SEB can diverge from the observations, but a suitable value for surface roughness can solve the issue (Fig. 6). When wind speed is high the SEB is especially sensitive to a change in roughness length, further highlighting the importance of wind related processes. At the ablation sites, the TI and ETI factors calibrated are highly variable, and depend on the season, site and year (Table S4). At ablation zone sites, pre-monsoon melt totals (Fig. 6, Table 4) are low and driven by SWnet. Factors are calibrated in order to reproduce observed ablation, even though daily mean Ta is not a good indicator of the pre-monsoon SEB (i.e. Results, Fig. 5, Table 3). Using individual factors for snow and ice directly includes the effect of changing surface α but results in increased factor variability. ETI shortwave radiation factors (SRF) are more stable from site to site and from year to year (Table S4). Ablation factors remain slightly variable, probably due to the presence of ablation processes such as surface sublimation, wind-driven surface erosion and snow redistribution. These processes do not relate well to Ta, contributing to changes of the TI or the ETI factors. TI factors for snow show low variability during the monsoon season, though this is likely due to low snow melt totals. Conversely, calibrated ice ablation factors range from 1 to 45 kg m−2 day−1 °C−1. At ablation zone sites, calibrated ice ablation factors are more variable and the calculated ablation is thus more sensitive to the parameter choices than during pre-monsoon.

Empirical ablation models need to accurately simulate the response of melt to large variations in SWnet. These changes primarily occur due to rapid and large changes in α caused by frequent precipitation events that deposit thin layers of snow on the ice surface. There is an additional inter-site and inter-year variability of TFTI and TFETI. Outgoing LW remains constant since Ts is always 0 during this melt-dominated period, so variations of cloudiness and LWout must be partly responsible for inter-year variability in calibrated ablation factors. Changes in the surface sublimation or wind-driven erosion and snow redistribution, which are not well related to Ta, could also contribute to ablation factor variability. The whole assessment is also strongly limited by the point-scale nature of our measurement set-up. In the accumulation zone of Mera Glacier, some very limited surface lowering was observed (Fig. 2) but was poorly captured by the TI, the ETI, and the ablation derived from SEB (ASEB, Fig. 6, see also Table S3 in Supplementary Material). The actual lowering can be explained only by surface compaction, densification or by wind-driven erosion and sublimation. Even though SWnet is the main input for ASEB it does not vary much since α remains high. Furthermore, small increases in SWnet only occur under clear skies, when LWnet is a large loss and w is stronger, which both favor low Ta. LWnet is more variable at the accumulation zone site than at lower elevations, and it drives changes in Ta.

At high elevations, the SEB relates poorly to Ta and SWnet, and furthermore represents only one part of mass loss. Therefore, the use of TI/ETI approaches to model ablation in this context is questionable. We conclude that for the presented sites, ablation models should (1) account for highly variable changes in albedo and cloudiness that affect melt rates, and (2) accommodate additional ablation processes such as wind-driven snow erosion and snow redistribution, as well as blowing snow sublimation. Improved and expanded observational studies at high elevation sites are thus necessary26. Models that can incorporate such processes (e.g. MESO-NH33,34 or COSMO-WRF35) should be explored with suitable test data. The correlation study conducted here shows that α has a strong control on the SEB. However, it is the change from snow to ice cover that exerts this strong control, since in the accumulation zone where surface remains snow covered, the correlation is much lower than in the ablation zone where surface transitions are more frequent. This highlights the additional importance of correctly characterizing the snow/rain temperature threshold. We note that α at the onset of monsoon is relevant for the total ablation over the season, and this is likely controlled by the amount of fallen spring snow.

Conclusions and Perspectives

In this study, meteorological controls on glacier ablation in High Mountain Asia (HMA) were studied using in-situ meteorological data collected at two glaciers in Nepal. We estimated the amount of mass loss due to surface melt and other ablation processes. The main variables controlling surface ablation were identified, and parameter values for a temperature-index (TI) and an enhanced temperature-index (ETI) model were provided. Modelling skill and parameters were assessed with respect to the observed meteorology at ablation zone and accumulation zone sites. Both models, calibrated for one season, perform well in the ablation zone during pre-monsoon and monsoon, since ablation is dominated by melt for those settings. A reliable estimate of surface α is nevertheless essential, since surface energy balance (SEB) changes are controlled by the changes in net shortwave radiation. Due to site and season-specific changes in the climatic drivers of ablation, the factors found for modeling ablation with a TI or an ETI are only transferable from one site to the other or from year to another year for what we define as the pre-monsoon period. During monsoon, surface melt (intense at times) dominates net ablation, though it has a site-dependent rate. TI and ETI calibrations are influenced by inter-annual changes in cloudiness. Also, site-specific changes in the snow-rain temperature point, in fresh-snow density, or in evaporation and sublimation rates, prevent the transferability of calibrated parameters to another site. During post-monsoon, melt remains limited and other ablation processes such as snow erosion by wind or sublimation are important but not statistically related to air temperature. This leads to the conclusion that a real alternative33,34,35 to TI and ETI should be developed for these cases. Existing models that are able to account for these processes32,33,35 are currently still computationally expensive. Surface lowering can be observed on some days at the high elevation site (Mera Glacier, 6352 m a.s.l.) during pre-monsoon and monsoon. This lowering is unlikely to be driven by melt, and is rather due to sublimation or other wind-driven processes, especially during pre-monsoon. Some ablation through surface melt may be occurring during monsoon but it remains weak and meltwater is typically refrozen in the snowpack below. The TI and ETI approaches cannot predict this weak high elevation monsoon melt, especially if calibrated in the ablation zone, as at high elevations air temperatures are below the necessary threshold for melt onset. The daily mean Ta relates poorly to melt totals as the SEB is positive mainly when there are clouds and when sublimation is reduced, due to lower wind speeds. The overall dependence of melt on α highlights the importance of collecting accurate α measurements and developing α models specific for high elevation sites. Accurate quantification of snowfall amounts and timing are consequently critical, and high elevation studies on the rain-snow threshold temperature are required.

Methods

Study sites

We used detailed on-glacier measurements between 2014 and 2017 on Yala Glacier (28.81 N, 85.84E, 1.61 km2) in the Langtang Valley and on Mera Glacier (27.7 N, 86.9E, 5.1 km2) in the Khumbu region. Measurements are made at accumulation and ablation zone sites at elevations ranging between 5350 and 6352 m a.s.l., (Fig. 1, Table 2). In addition, we use measurements from two AWS fixed on bedrock near the terminus of each glacier (Fig. 1, Table 2), and data from nearby pluviometers in the Khumbu valley. Mera Glacier, located approximately 30 km south of Mt. Everest, straddles the Hinku Valley and Hunku Valley in Central Nepal. There are two branches of the same glacier in the ablation zone, known as Naulek Glacier and Mera Glacier. Mera Glacier ranges in elevation from 4940 to 6420 m a.s.l.18, faces north in the accumulation area and south-east in the ablation area of Naulek branch. The monitoring program on Mera Glacier was initiated in 200718,28. Yala Glacier is a small south-west facing glacier located approximately 135 km west of Mera Glacier in Langtang Valley, Nepal. Yala Glacier ranges in elevation from 5128 to 5749 m a.s.l.23,36.

Data

On-glacier meteorological measurements were obtained from ‘floating’ stations (Yala Glacier, May to October 2014) and stations drilled into the ice surface (Yala Glacier, May 2016 to October 2017). At Mera Glacier, 5380 m a.s.l., ablation zone site, an AWS operated continuously between November 2012 and November 2016, while the AWS on Mera Glacier, 6352 m a.s.l. (near summit), was operational intermittently between November 2013 and July 2016 due to difficult meteorological conditions from 26 August 2014. Continuous meteorological data from Mera Glacier, 6352 m a.s.l., are only available from November 2013 to 26 August 2014, from December 2014 to March 2015, from April 2015 until August 2015 then December 2015 to 2 August 2016 (Table 1). The Mera Glacier, 5380 m a.s.l., ablation zone station, was drilled in ice until March 2014, and after was on a free standing tripod. Mera Glacier, 6352 m a.s.l., station was drilled in ice. Additionally, at both glacier sites in the Langtang and Everest areas, an AWS was permanently operated nearby the glacier tongue on rocky and ice-free surfaces, at 5350 m a.s.l. on Mera La for the Naulek-Mera site and at 5090 m a.s.l at Yala Glacier base camp. All stations measure the 4 components of the surface radiative balance, air temperature and humidity, wind speed and direction, and surface height changes resulting a priori from ice or from snow deposition or ablation processes. Precipitation is recorded at Yala Glacier base camp AWS. Precipitation was also measured in the Khumbu region at the Pyramide site and at Pheriche (Table 1, Fig. 1). Atmospheric pressure was measured intermittently at the different AWS locations. Details of the settings and data availability can be found in Table 2.

Surface energy balance and statistics

We defined the hourly surface energy balance (SEB) in W m−2 as:

$${\rm{S}}{\rm{E}}{\rm{B}}=(S{W}_{inc}+S{W}_{out}+L{W}_{inc}+L{W}_{out}+H+LE)$$
(1)

Where SWinc, SWout, LWinc and LWout are the incoming and outgoing shortwave and longwave radiation, H and LE are the turbulent fluxes of sensible and latent heat fluxes, respectively. The individual SEB components were measured or calculated from the AWS data. The radiative fluxes were given by the radiometers. The turbulent fluxes H and LE were derived using the so-called bulk-aerodynamic method37, and using the glacier AWS measurements: barometric pressure (Pa), wind speed (w), air temperature (Ta) and relative humidity (RH) and snow/ice surface temperature Ts calculated from the measured surface longwave emission, assuming snow and ice emissivity (ε) was equal to 0.99. The aerodynamic (z0), thermal (z0t) and humidity (z0q) roughness lengths were set constant to 0.01 m, 0.001 m and 0.001 m, respectively. We used stability corrections as in Litt et al.37. The values of the aerodynamic roughness length can be seen as high but are based on Stigter et al.26, who calibrated them on Yala Glacier so that bulk fluxes fit fluxes as evaluated from an eddy-covariance tower.

In order to provide insights on the meteorological processes governing SEB, we calculated the Pearson correlation coefficient (R) between (i) each daily mean component of the SEB and the daily mean SEB, (ii) of the temperature with the daily mean SEB and (iii) with all the daily mean SEB components. We then (iv) calculate the contribution of each SEB term x to the correlation between the SEB and the daily mean Ta and (v) to the correlation between the SEB and the shortwave radiation, as:

$${\rm{R}}{({T}_{{\rm{a}}},{\rm{SEB}})}_{x}=\frac{{\sigma }_{x}}{{\sigma }_{SEB}}{\rm{R}}({T}_{a},x)$$
(2)
$${\rm{R}}{(S{W}_{net},{\rm{SEB}})}_{x}=\frac{{\sigma }_{x}}{{\sigma }_{SEB}}{\rm{R}}(S{W}_{{\rm{net}}},x)$$
(3)

These partial correlations indicate the degree of variability of variable x, relating to Ta (Eq. 2) or to the SWnet (Eq. 3), and also to the net SEB. If (2) or (3) are high for a given x and R(Ta, SEB) is high, this indicated that the x variable explained well the SEB variability, in relation to Ta.

Ablation models

Since we did not measure the fluxes below the surface, we used hourly surface temperature as a proxy for melt occurrence. When Ts reached 0 °C once after the winter, hourly negative SEB values were accumulated as a “cold storage” term. This cold storage had to be compensated for and as soon as the SEB was positive and surface temperature was 0 °C we converted the energy to hourly melt. The sum of the melt derived from the SEB plus the ablation due to sublimation and evaporation derived from LE was called ASEB in kg m−2. This was validated against the observed ablation (Aobs). It was obtained by converting the measured daily surface height changes to water equivalent, using a density (ρ, kg m−3) parameterization based on the observed α, inspired by that used in38 and tuned to the observations such as:

$$\alpha < 0.2,\rho =850\,{\rm{kg}}\,{{\rm{m}}}^{-3}({\rm{ice}})$$
(4)
$$\alpha > 0.2,\rho =680\to 5\,{\rm{kg}}\,{{\rm{m}}}^{-3}({\rm{firn}})$$
(5)

Temperature index (TI) and enhanced temperature index (ETI) Models

To account for specificities of the site and season, we calibrated a unique set of the TI and ETI parameters. We discuss its ability to model observed ablation at other sites, comparing the calibrated parameters to each other (Table S4). For the TI we have, if the daily mean of Ta is larger than TTI:

$${A}_{{\rm{TI}}}=T{F}_{{\rm{TI}}}(\overline{{T}_{{\rm{a}}}}-{T}_{{\rm{TI}}})$$
(6)

The overbar indicates daily mean. TFTI kg m−2 °C −1 day−1 is given a different value when α is higher than 0.4 (snow-melt factor) or lower than 0.4 (ice-melt factor). For the ETI model we have, if the daily mean Ta is larger than TETI:

$${A}_{{\rm{ETI}}}=T{F}_{{\rm{ETI}}}(\overline{{T}_{{\rm{a}}}}-{T}_{{\rm{ETI}}})+SRF.S{W}_{{\rm{inc}}}(1-\alpha )$$
(7)

The TFETI is set the same for all surface states, and AETI stands for the daily ablation (m we) derived from the SRF (m we W−1 m2 day−1) and the daily mean SWinc. We calibrate the two parameters TFTI for snow and for ice, the TFETI, and the SFR by generating an ablation time series for each possible combination of these parameters and choose the one to minimize the factor f as defined below:

$$f={\rm{\Sigma }}|\frac{{A}_{{\rm{cum}}}-{A}_{{\rm{obs}}}}{{A}_{{\rm{obs}}}}|$$
(8)

Where Acum is the cumulated ATI or AETI, and Aobs is derived with the sounding height ranger and the density parameterization (Eqs 4 and 5), at a daily time scale. The TF values are between 0 and 100 kg m−2 °C−1 and SRF between 0 and 1 kg m−2 W−1 m2 day−1. We considered the calibration failed if, to minimize f, one of the parameters must take its maximum value, or if all parameters must be 0.